# OperatorSchmidtDecomposition

 Other toolboxes required OperatorSchmidtDecomposition Computes the operator Schmidt decomposition of a bipartite operator none FilterNormalFormIsProductOperatorOperatorSchmidtRankSchmidtDecompositionTensorSum Entanglement and separability

OperatorSchmidtDecomposition is a function that computes the operator Schmidt decomposition of a bipartite operator. The user may specify how many terms in the operator Schmidt decomposition they wish to be returned.

## Syntax

• S = OperatorSchmidtDecomposition(X)
• S = OperatorSchmidtDecomposition(X,DIM)
• S = OperatorSchmidtDecomposition(X,DIM,K)
• [S,U,V] = OperatorSchmidtDecomposition(X,DIM,K)

## Argument descriptions

### Input arguments

• X: A bipartite operator to have its operator Schmidt decomposition computed.
• DIM (optional, by default has both subsystems of equal dimension): A specification of the dimensions of the subsystems that X lives on. DIM can be provided in one of three ways:
• If DIM is a scalar, it is assumed that the first subsystem has dimension DIM and the second subsystem has dimension length(X)/DIM.
• If $X \in M_{n_1} \otimes M_{n_2}$ then DIM should be a row vector containing the dimensions (i.e., DIM = [n_1, n_2]).
• If the subsystems aren't square (i.e., $X \in M_{m_1, n_1} \otimes M_{m_2, n_2}$) then DIM should be a matrix with two rows. The first row of DIM should contain the row dimensions of the subsystems (i.e., the mi's) and its second row should contain the column dimensions (i.e., the ni's). In other words, you should set DIM = [m_1, m_2; n_1, n_2].
• K (optional, default 0): A flag that determines how many terms in the operator Schmidt decomposition should be computed. If K = 0 then all terms with non-zero operator Schmidt coefficients are computed. If K = -1 then all terms (including zero terms) are computed. If K > 0 then the K terms with largest operator Schmidt coefficients are computed.

### Output arguments

• S: A vector containing the operator Schmidt coefficients of X.
• U (optional): A cell containing the left Schmidt operators of X.
• V (optional): A cell containing the right Schmidt operators of X.

## Examples

### A random example

The following code computes the operator Schmidt decomposition of a random density matrix in $M_2 \otimes M_4$ and verifies that it is indeed a valid tensor decomposition of that density matrix:

>> rho = RandomDensityMatrix(8);
>> [s,U,V] = OperatorSchmidtDecomposition(rho,[2,4])

s =

0.3986
0.2310
0.1637
0.1118

U =

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

V =

[4x4 double]    [4x4 double]    [4x4 double]    [4x4 double]

>> norm(s(1)*kron(U{1},V{1}) + s(2)*kron(U{2},V{2}) + s(3)*kron(U{3},V{3}) + s(4)*kron(U{4},V{4}) - rho)

ans =

1.8449e-016

As an alternative (and easier) way to verify that we have a valid tensor decomposition of rho, we can use the TensorSum function:

>> norm(TensorSum(s,U,V) - rho)

ans =

1.8449e-016

## Source code

Click on "expand" to the right to view the MATLAB source code for this function.

1. %%  OPERATORSCHMIDTDECOMPOSITION   Computes the operator Schmidt decomposition of a bipartite operator
2. %   This function has one required argument:
3. %     X: a bipartite operator
4. %
5. %   S = OperatorSchmidtDecomposition(X) is a vector containing the non-zero
6. %   operator Schmidt coefficients of the bipartite operator X, where the
7. %   two subsystems are each of size sqrt(length(X)) and X is assumed to be
8. %   square.
9. %
10. %   This function has two optional input arguments:
11. %     DIM (default has both subsystems of equal size)
12. %     K (default 0)
13. %
14. %   [S,U,V] = OperatorSchmidtDecomposition(X,DIM,K) gives the operator
15. %   Schmidt coefficients S of the operator X and the corresponding left and
16. %   right tensor operators in the cells U and V. DIM is a 1x2 vector containing
17. %   the dimensions of the subsystems that X lives on. If X is non-square,
18. %   different row and column dimensions can be specified by putting the row
19. %   dimensions in the first row of DIM and the column dimensions in the
20. %   second row of DIM. K is a flag that determines how many terms in the
21. %   operator Schmidt decomposition should be computed. If K = 0 then all
22. %   terms with non-zero operator Schmidt coefficients are computed. If
23. %   K = -1 then all terms (including zero terms) are computed. If K > 0
24. %   then the K terms with largest operator Schmidt coefficients are
25. %   computed.
26. %
27. %   URL: http://www.qetlab.com/OperatorSchmidtDecomposition
28. 
29. %   requires: opt_args.m, PermuteSystems.m, SchmidtDecomposition.m, Swap.m
30. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
31. %   package: QETLAB
32. %   last updated: November 12, 2014
33. 
34. function [s,U,V] = OperatorSchmidtDecomposition(X,varargin)
35. 
36. dX = size(X);
37. sdX = round(sqrt(dX));
38. 
39. % set optional argument defaults: dim=sqrt(length(X)), k=0
40. [dim,k] = opt_args({ [sdX(1) sdX(1);sdX(2) sdX(2)], 0 },varargin{:});
41. 
42. % allow the user to enter a single number for dim
43. if(length(dim) == 1)
44.     dim = [dim,dX(1)/dim];
45.     if abs(dim(2) - round(dim(2))) >= 2*dX(1)*eps
46.         error('OperatorSchmidtDecomposition:InvalidDim','If DIM is a scalar, X must be square and DIM must evenly divide length(X); please provide the DIM array containing the dimensions of the subsystems.');
47.     end
48.     dim(2) = round(dim(2));
49. end
50. 
51. % allow the user to enter a vector for dim if X is square
52. if(min(size(dim)) == 1)
53.     dim = dim(:)'; % force dim to be a row vector
54.     dim = [dim;dim];
55. end
56. 
57. % The operator Schmidt decomposition is just the Schmidt decomposition of a
58. % related vector obtained by moving matrix elements around.
59. [s,u,v] = SchmidtDecomposition(Swap(reshape(X,prod(prod(dim)),1),[2,3],[dim(2,:),dim(1,:)]),prod(dim),k);
60. 
61. % Now reshape things into the proper output format.
62. sr = length(s);
63. U = mat2cell(reshape(u,dim(1,1),dim(2,1)*sr),dim(1,1),dim(2,1)*ones(1,sr));
64. V = mat2cell(reshape(v,dim(1,2),dim(2,2)*sr),dim(1,2),dim(2,2)*ones(1,sr));