# 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? Sep 18, 2020 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. Sep 19, 2020 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 Sep 21, 2020 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). Sep 21, 2020 at 22:28
• Hm, thanks, indeed that fixes it. Not that I understand how it is possible... Sep 22, 2020 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. Sep 19, 2020 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. Sep 20, 2020 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]? Sep 20, 2020 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. Sep 20, 2020 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. Sep 21, 2020 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. Sep 19, 2020 at 11:44