Логарифм сингулярных матриц

Question

Я хотел рассчитать следующее количество:

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

как в относительной энтропии . Вот,$\rho_1, \rho_2$ положительные полуопределенные матрицы с неортогональным носителем (чтобы не расходиться - такие термины, как $0 \log[0]$ должно быть сведено к нулю, что является стандартным предположением) и $Tr$ это след.

Есть аналогичный вопрос SE , данная функция

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

который должен иметь дело с матричным логарифмом, ведет себя, однако, странно.

Например, я предполагаю, что $\rho_1 =\rho_2 =\{\{0.33,0,0\},\{0,0,0\},\{0,0,0.66\}\}$. Количество$X$ тогда должно быть

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

Однако использование MatriLogSafe в определении дает другой результат:

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

Проблема в том, что MatrixLogSafe иногда «переключает собственные векторы»,

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}}

(так $\log[0.33]= -1.10866$ и $\log[0.66]=-0.415515$, но вывод должен быть {{-1.10866, 0, 0.}, {0., 0., 0.}, {0., 0., - 0.415515}}).

(Почему-то я думаю, что проблема в том, что я использую числовые значения, но я хочу, чтобы функция работала как для числовых, так и для «точных» (?) Чисел)

Как это исправить?

---

Я рассмотрел ответ Карла Уолля, но все еще что-то не работает. В частности, рассмотрим две матрицы,$\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}}

и $\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])}}

которые могут быть уродливыми, но они оба эрмитовские и положительно-полуопределенные.

Я как бы сложен $X$ ($X=densityTrace[\rho_1,\rho_2]=-0.0019613 + 0.393667 I$), чего не может быть!

Однако, если я сначала рассчитываю числовые матрицы, я получаю реальный результат $densityTrace[N[\rho_1],N[\rho_2]]=-0.0432473$ (?!)

Я немного изменил определение 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, λ}]

]

(добавляет Чоп), чтобы избавиться от очень маленьких воображаемых «отходов», это хорошая идея?

Answer 1

Может быть, вы можете использовать Eigensystemдля преобразования своей сингулярной матрицы$\rho_2$ в $S.d.S^{-1}$, а затем использовать циклическое свойство Tr. То есть:

$$ \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} $$

Код, реализующий эту идею:

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

Answer 2

Наверное, проще обнаружить сингулярность раньше:

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

Я могу ошибаться относительно того, какое именно состояние вам нужно. Во всяком случае, я не ожидаю, что это принесет MatrixFunction[Piecewise[...]]вам много пользы. Я даже не уверен, что это условие # > 0должно делать внутри MatrixFunction.