The operator-Schmidt decomposition does not work on Hermitian operators

[denisrosset]

Just realized that the operator-Schmidt decomposition function OperatorSchmidtDecomposition works on the basis of arbitrary operators, but Hermitian ones.

This runs contrary to the definition e.g. in http://www.njohnston.ca/2014/06/what-the-operator-schmidt-decomposition-tells-us-about-entanglement/

One way to perform this decomposition is to first decompose the Hermitian operator on a basis of (local) Hermitian matrices such that

rho = \sum_{ij} c_{ij} A_i (x) B_j

where the c_{ij} are real, perform the SVD on c_{ij}, and then expand back the orthonormal Hermitian operators from the result [U S V] = svd(c). This is how I'm doing it for myself, but I'd like to submit a PR for QETLAB.

Is there a better way?

[nathanieljohnston]

Can you provide an explicit example of it not working on Hermitian operators? Then (random) examples that I try of Hermitian matrices return Hermitian terms in the operator Schmidt decomposition, but I can see how it might mess up if some singular values are repeated (in which case it indeed needs fixing). Thanks!

[denisrosset]

Yep! Thanks for your quick reply.

Note: I am curious about implementing this with reshape/permute tricks, as I could not understand exactly what the code is doing.

>> [c A B] = OperatorSchmidtDecomposition([1 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 1], [2 2])
[c A B] = OperatorSchmidtDecomposition([1 0 0 1; 0 0 0 0; 0 0 0 0; 1 0 0 1], [2 2])

c =

     1
     1
     1
     1


A = 

    [2x2 double]    [2x2 double]    [2x2 double]    [2x2 double]


B = 

    [2x2 double]    [2x2 double]    [2x2 double]    [2x2 double]

>> A{2}
A{2}

ans =

     0     0
     1     0

[nathanieljohnston]

Thanks! Yeah, as I suspected the problem here is repeated singular values in the SVD. That'll be a bit of a pain to get around in a manner as efficient as I'd like, so I'll have to think a bit more about what the "best" way to fix this would be (I really don't want to do 2 SVDs.

The reshape tricks are basically just writing the matrix in a particular basis so that computing that SVD of the resulting reshaped matrix is the operator-schmidt decomp of the matrix we actually care about. If memory serves, we're just converting the basis matrix |i><j| \otimes |k><\ell| to |i><k| \otimes |j><\ell|. To do this just with Hermitian matrices we probably want to do a similar trick but with some usage of triu thrown in so that we only look at the upper-triangular entries and write those in a Hermitian basis.

[denisrosset]

Thanks! indeed, that can be a bit of a pain.