# Logarithm of singular matrices

I wanted to calculate following quantity:

$$X=Tr[\rho_1 \log[\rho_2]]$$,

as in the relative entropy. Here, $$\rho_1, \rho_2$$ are positive semidefinite matrices with non-orthogonal support (so that the thing does not diverge - terms like $$0 \log[0]$$ should be taken to zero, which is standard assumption) and $$Tr$$ is trace.

There is a similar SE question, the given function

MatrixLogSafe[x_] := MatrixFunction[Piecewise[{{Log[#1], #1 > 0}}] &, x]

which should deal with the matrix logarithm, it behaves, however, strange.

For e.x. I assume that $$\rho_1 =\rho_2 =\{\{0.33,0,0\},\{0,0,0\},\{0,0,0.66\}\}$$. The quantity $$X$$ should then be

$$X=0.33 \log[0.33] + 0.66 \log[0.66]= -0.640099$$.

However, using MatriLogSafe in the definition gives different output:

In[402]:= Tr[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}.MatrixLogSafe[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}]]

Out[402]= -0.731717

The problem is, that MatrixLogSafe sometimes "switch the eigenvectors",

In[403]:= MatrixLogSafe[{{0.33, 0, 0}, {0, 0, 0}, {0, 0, 0.66}}]

Out[403]= {{0., 0., 0.}, {0., -0.415515, 0.}, {0., 0., -1.10866}}

(so $$\log[0.33]= -1.10866$$ and $$\log[0.66]=-0.415515$$, but the output should be { { -1.10866, 0, 0.}, {0., 0., 0.}, {0., 0.,-0.415515}}).

(Somehow I think the problem is that the I use numerical values, but I want that the function works for both numerical and "exact" (?) numbers)

How one can fix it?

I have consider the answer given by Carl Woll, however still something is not working. In particular, consider two matrices, $$\rho_1$$:

{{1/4, 1/4 E^(-((I \[Pi])/10)), 1/4 E^(-((I \[Pi])/10)), 1/
4}, {1/4 E^((I \[Pi])/10), 1/4, 1/4,
1/4 E^((I \[Pi])/10)}, {1/4 E^((I \[Pi])/10), 1/4, 1/4,
1/4 E^((I \[Pi])/10)}, {1/4, 1/4 E^(-((I \[Pi])/10)),
1/4 E^(-((I \[Pi])/10)), 1/4}}

and $$\rho_2$$:

{{1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]),
1/4 E^(-(1/2) - (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)),
1/4 E^(-(1/2) - (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)),
1/(4 Sqrt[
E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))^2}, {1/4 E^(-(1/2) + (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)), 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]),
1/(4 Sqrt[
E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))^2,
1/4 E^(-(1/2) + (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))}, {1/4 E^(-(1/2) + (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)),
1/(4 Sqrt[
E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))^2, 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E]),
1/4 E^(-(1/2) + (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))}, {1/(
4 Sqrt[E]) + (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10))^2,
1/4 E^(-(1/2) - (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)),
1/4 E^(-(1/2) - (I \[Pi])/10) +
1/2 (1 - 1/Sqrt[E]) (1/4 E^(-((I \[Pi])/10)) +
1/4 E^((I \[Pi])/10)), 1/4 (1 - 1/Sqrt[E]) + 1/(4 Sqrt[E])}}

which are maybe ugly, but they both are Hermitian and positive-semi-definite.

I get however complex $$X$$ ($$X=densityTrace[\rho_1,\rho_2]=-0.0019613 + 0.393667 I$$), which cannot be!

However, if I calculate firstly the numerical matrices I get real output $$densityTrace[N[\rho_1],N[\rho_2]]=-0.0432473$$ (?!)

I have slightly changed the definition of densityTrace

densityTrace[a_, b_] := Module[{λ, S, d},
{λ, S} = Eigensystem[b];
S = Transpose[S];
d = Diagonal[Inverse[S]. a. S];
Total @ MapThread[If[Chop[#1]==0,0,Chop[#1] Log[Chop[#2]]]&, {d, λ}]

]

(adding Chop), to get rid of some very small imaginary "waste", is it a good idea?

• Why not find the eigenvalues of $\rho_1$ instead, and drop zero eigenvalues? – Carl Woll Sep 18 '20 at 15:58
• Can you be more specific? In this case when $\rho_1=\rho_2$ it would works, but not for $\rho_1\neq\rho_2$, if I get you right. – Agnieszka Sep 19 '20 at 11:54
• sounds like you are taking log of of some negative numbers. What's inside log(d)? Also your example is hard to reproduce, better give example in Mathematica code in addition to latex – Yaroslav Bulatov Sep 21 '20 at 15:20
• For complicated symbolic matrices, it is better to use extended precision arithmetic instead of machine number arithmetic. Compare densityTrace[N[r1, 50], N[r2, 50]] with N[densityTrace[r1, r2], 50] and you will see that they basically agree (note that you shouldn't be using Chop in the definition of densityTrace). – Carl Woll Sep 21 '20 at 22:28
• Hm, thanks, indeed that fixes it. Not that I understand how it is possible... – Agnieszka Sep 22 '20 at 10:17

Maybe you can use Eigensystem to transform your singular matrix $$\rho_2$$ into $$S.d.S^{-1}$$, and then use the cyclical property of Tr . That is:

\begin{align} \text{Tr}{\left(\rho_1 . \log(\rho_2) \right)}&= \text{Tr}{\left(\rho_1 . \log(S . d . S^{-1})\right)} \\ &= \text{Tr}{\left(\rho_1 . S . \log(d) . S^{-1} \right)}\\ &= \text{Tr}{\left(S^{-1} . \rho_1 . S . \log(d)\right)}\\ &= \text{Diagonal}{\left( S^{-1} . \rho_1 . S\right) . \log (d)} \end{align}

Some code that implements this idea:

densityTrace[a_, b_] := Module[{λ, S, d},
{λ, S} = Eigensystem[b];
S = Transpose[S];
d = Diagonal[Inverse[S]. a. S];
Total @ MapThread[If[#1==0,0,#1 Log[#2]]&, {d, λ}]
]

m1 = m2 = DiagonalMatrix[{.33, 0, .66}];
densityTrace[m1, m2]

-0.640099

• Ok, thanks, I think this will work. – Agnieszka Sep 19 '20 at 17:45
• Nice, but Total@MapThread[ If[#1 == 0., 0, #1 Log[#2]] &, {{0.66, 0.33, 0.}, {0.66, 0.33, 0.}}]` would suffice for this kind of question. And is faster, Theory is theory and that is good. – Steffen Jaeschke Sep 20 '20 at 9:57
• @SteffenJaeschke How would your faster version work if the matrices were With[{s = RandomReal[1, {3, 3}]}, Inverse[s].DiagonalMatrix[{.33, 0, .66}].s]? – Carl Woll Sep 20 '20 at 16:35
• Excuse. I just considered the given example strictly. A density composed of three real, positive number of scalar type put together and possibly both the very same. Then my suggestion is use a type check m_?MatrixQ or v_?VectorQ and cast type for professionality. Thanks. – Steffen Jaeschke Sep 20 '20 at 17:16
• @Agnieszka You can't possibly expect anyone to type those matrices into Mathematica! You need to include a copyable version in your question. – Carl Woll Sep 21 '20 at 15:47

It's probably easier to detect the singularity earlier:

MatrixLogSafe[x_?SquareMatrixQ] /; Det[x] > 0 := MatrixLog[x];
MatrixLogSafe[x_?SquareMatrixQ] := ConstantArray[0, Dimensions[x]];

I could be wrong about the exact condition you need. At any rate, I don't expect MatrixFunction[Piecewise[...]] to do you much good. I'm not even sure what the condition # > 0 is supposed to do inside of MatrixFunction.

• Hm, this does not work, if I get it right for a singular matrix it sets it simply to zero (I mean the whole matrix), so that in my particular example of X the output is 0. – Agnieszka Sep 19 '20 at 11:44