# How to calculate the trace of a product of matrices taking advantage of its properties?

Hello I'm trying to compute the trace of the product of matrices like this

$$$$\begin{pmatrix} a+\imath b &c &0 \\ c& d &c \\ 0& c &f+\imath g \end{pmatrix}^{-1} \begin{pmatrix} -2b &0 &0 \\ 0 & 0 &0 \\ 0& 0 &0 \end{pmatrix} \begin{pmatrix} a+\imath b &c &0 \\ c& d &c \\ 0& c &f+\imath g \end{pmatrix}^{*-1} \begin{pmatrix} 0 &0 &0 \\ 0 & 0 &0 \\ 0& 0 &-2g \end{pmatrix}$$$$

where $*$ is the comlpex transpose and all the components are real.

This situation is very easy to compute it but if I want to do it for a 50x50 matrix (the first is a tridiagonal matrix with $d$ in the main diagonal except the $(1,1)$ and $(N,N)$ elements, $c$ in the other two diagonals,the second and the forth have non zero elements only in the $(1,1)$ and $(N,N)$ positions) it's almost take a lot of time.

Taking advantage of the symmetry and that at the end only the last column is non zero (and therefore only the element $(N,N)$ will contribute) is there a way to compute it faster?

Thanks

Edit: The code I'm using for a 10x10 matrix is the following

    M = {{a + I b, c, 0, 0, 0, 0, 0, 0, 0, 0}, {c, d, c, 0, 0, 0, 0, 0, 0,
0}, {0, c, d, c, 0, 0, 0, 0, 0, 0}, {0, 0, c, d, c, 0, 0, 0, 0,
0}, {0, 0, 0, c, d, c, 0, 0, 0, 0}, {0, 0, 0, 0, c, d, c, 0, 0,
0}, {0, 0, 0, 0, 0, c, d, c, 0, 0}, {0, 0, 0, 0, 0, 0, c, d, c,
0}, {0, 0, 0, 0, 0, 0, 0, c, d, c}, {0, 0, 0, 0, 0, 0, 0, 0, c,
f + I g}};

G1 = {{-2 b, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0,
0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0}};
G2 = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0,
0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0,
0, -2 g}};

Inverse[M].G1.Inverse[ConjugateTranspose[M]].G2 //ComplexExpand


To make tridiagonal matrices I'm using the comand

    SparseArray[{Band[{1, 2}] -> c, Band[{1, 1}] -> d, Band[{2, 1}] -> c},10] // Normal


and then recopying it, change only the two elements I need to change and then I operate

• Do you want a symbolic answer in terms of $a,b,c,d,f,g$, or are you happy with a function that takes numerical values for them and returns the trace? The two will obviously perform rather differently. – Emilio Pisanty Aug 12 '16 at 14:31
• @EmilioPisanty symbolic answer because each of this terms are functions of other variables but always this varables appear on the same way – Daniel Aug 12 '16 at 14:35
• Can you enter your matrices in copy-and-paste-able Mathematica code, properly formatted in code blocks? It will make it much easier for use to help you if we can just copy and paste into our own copies of Mathematica rather than have to type it in ourselves. – march Aug 12 '16 at 15:35
• @march I edited the post with the code I'm using for a 10x10 matrix – Daniel Aug 12 '16 at 15:45
• I'll think about this, but it might be worth writing this using the Hadamard product (i.e. element-wise multiplication of matrices) in this form. – march Aug 12 '16 at 15:55

The following exploits sparse matrices, and uses LinearSolve to avoid computing the inverse which is not sparse. I suspect that further optimisation would be possible, but it is certainly practical for problems of the scale you describe. (I have assumed that you requre a numeric, rather than a symbolic solution)

test[n_] :=  Module[{c = 1.3, b = 0.2, a = 4.1, f = 2.2, g = -0.3, d = 0.5, M, G1, G2},
M = SparseArray[{Band[{1, 2}] -> c, Band[{1, 1}] -> d, Band[{2, 1}] -> c}, n];
M[[1, 1]] = a + I b;
M[[-1, -1]] = f + I g;
G1 = SparseArray[{1, 1} -> -2 b, {n, n}];
G2 = SparseArray[{-1, -1} -> -2 g, {n, n}];
Tr[LinearSolve[M, G1].LinearSolve[ConjugateTranspose[M], G2]]]


In my tests, this agrees with the formula that you gave. I have

AbsoluteTiming[test[5000]]
(* {9.70897, -0.0233494 + 7.80626*10^-17 I} *)