OperatorSchmidtDecomposition

From QETLAB
Jump to: navigation, search
OperatorSchmidtDecomposition
Computes the operator Schmidt decomposition of a bipartite operator

Other toolboxes required none
Related functions FilterNormalForm
IsProductOperator
OperatorSchmidtRank
SchmidtDecomposition
TensorSum
Function category 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));