Is there an algorithm or graph theory that allows me to not need to store an intermediate matrix when calculating AT*Y1*A + BT*Y2*B?

Question

I have a system of conductors for which there are two dense matrices of the (complex) mutual admittances, $Y_A$ and $Y_B$, which are symmetric. Then, an equivalent nodal admittance matrix $Y_N$ is calculated by the following:

$$ Y_N = A^T \cdot Y_A \cdot A + B^T \cdot Y_B \cdot B $$

where $A$ is a directed incidence matrix (directed graph) for which each element $a_{ik}$ is given by

$$ a_{ik} = \begin{cases} 1 \quad \text{if node $k$ is the start point of conductor $i$} \\ -1 \quad\text{if node $k$ is the end point of conductor $i$} \\ 0 \quad\text{otherwise} \\ \end{cases} $$

and $B$ is an undirected incidence matrix for which the elements $b_{ik}$ are given by

$$ b_{ik} = \frac{|a_{ik}|}{2} = \begin{cases} 1 / 2 \quad \text{if node $k$ is connected to conductor $i$} \\ 0 \quad\text{otherwise} \\ \end{cases} $$

Each coductor has exactly one start point and exactly one end point, which makes $A$ and $B$ sparse, with only two non-zero entries per line. But I store them as dense matrices so I can pass them to BLAS and LAPACK routines.

The system have $m$ conductors and $n$ nodes. $Y_A$ and $Y_B$ are $m \times m$ and $A$ and $B$ are $m \times n$.

---

I am using BLAS to do those calculations by calling CSYMM and CGEMM two times each. The following pseudo-code summarizes the steps.

calculate A and B, storing them as dense complex matrices
for i = 1:N
    calculate YA and YB
    C := YA * A + 0 * C  # CSYMM
    YN := A^T * C + 0 * YN  # CGEMM
    C := YB * B + 0 * C  # CSYMM
    YN := B^T * YB + 1 * YN  # CGEMM
end

The system I want to simulate have about 50000 conductors and nodes. I am storing 6 dense matrices, $A, B, C, Y_N, Y_A, Y_B$, with 2.5 billions elements of complex float each. That requires 120 GB of memory and I would like to reduce that.

At first I thought about storing $A$ and $B$ as dense bit matrices (inspired by this question on StackOverflow), get the source code of cgemm and csymm at the netlib site and rewrite a version of them to use the bit matrices, doing the bit to complex type conversion inside the loop. That saves me almost 40 GB of memory and seems like a good solution, but I am wondering: is there an algorithm or graph theory that would allow me to not need to store the intermediate complex matrix $C$, saving me another 20 GB of memory?

Answer 1

BLAS may not have a function to compute what you are asking for, but the product $$ Y_N = A^TY_AA + B^T Y_B B $$ means that the entries $(Y_N)_{ij}$ are defined by $$ (Y_N)_{ij} = \sum_{k,l} (A^T)_{ik}(Y_A)_{kl}A_{lj} + \sum_{k,l} (B^T)_{ik}(Y_B)_{kl}B_{lj} $$ which is equivalent to $$ (Y_N)_{ij} = \sum_{k,l} \left\{ a_{ki}(Y_A)_{kl}a_{lj} + b_{ki}(Y_B)_{kl}b_{lj} \right\}. $$ Given your formula for how $A$ and $B$ are related, this is in turn equal to $$ (Y_N)_{ij} = \sum_{k,l} \left\{ a_{ki}(Y_A)_{kl}a_{lj} + \frac 14 |a_{ki}|(Y_B)_{kl}|a_{lj}| \right\}. $$

You can implement this operation using four nested loops over $i,j,k,l$, and if you make the loops over $i,j$ the outer ones, you don't need any intermediate storage -- in fact, you don't need to store the $B$ matrix either based on the formula. You really don't need to store the $A$ matrix either; all you need to know is whether $a_{ij}$ is zero, plus one, or minus one. That can be stored with far fewer bits than as a floating point number.

This quadruple loop may, if you do it naively, not be as fast as what BLAS does. But it saves memory, and that was after all what you were after. On the other hand, reducing memory requirements on modern processors often translates into faster computations, so it may go either way.

Answer 2

I agree with the comments that your best bet is likely to use a sparse multiplication so that you don't need to store A and B as dense.

If you want something easy that uses high level lapack/blas routines, consider the following algorithm:

1. Calculate cholesky decompositions $Y_a = L_aL_a^H$ and $Y_b = L_bL_b^H$ (You can use $Y_n$ for workspace).
2. $A := L_a^HA$ and $B := L_b^HB$ (ctrmm)
3. $Y_n := A^HA$ (cherk)
4. $Y_n := B^HB + Y_n$ (cherk)

Some downsides to this approach:

• Assumes matrices are positive definite, if not, you need to use LDLT instead
• Likely introduces extra rounding error
• Cholesky decomposition may be prohibitively expensive to calculate (doesn't scale as well as gemm)

Some upsides to this approach:

• Less memory, you can even store the triangular parts of $Y_a$ and $Y_b$ in the same $m \times m$ matrix if you store the diagonal separately.
• Likely faster than dense for loop approach. I think the for loop will have many cache misses. (And a quadruple for loop would be O(n^4), pretty expensive for matrices that large)