UPB

From QETLAB
Jump to: navigation, search
UPB
Generates an unextendible product basis

Other toolboxes required none
Related functions IsUPB
MinUPBSize
UPBSepDistinguishable
Function category Unextendible product bases

UPB is a function that generates an unextendible product basis (UPB). The user may either request a specific UPB from the literature such as 'Tiles' or 'Pyramid'[1], or they may request a minimal UPB of specified dimensions.

Syntax

  • U = UPB(NAME)
  • U = UPB(NAME,OPT_PAR)
  • [U,V,W,...] = UPB(NAME,OPT_PAR)
  • U = UPB(DIM)
  • U = UPB(DIM,VERBOSE)
  • [U,V,W,...] = UPB(DIM,VERBOSE)

Argument descriptions

Input arguments

Important: Do not specify both NAME and DIM: just one or the other!

  • NAME: The name of a UPB that is found in the literature. Accepted values are:
    • 'GenShifts': A minimal UPB in $(\mathbb{C}^2)^{\otimes p}$ (only valid when p ≥ 3 is odd) constructed in [2]. Note that OPT_PAR must be the number of parties (i.e., the integer p) in this case.
    • 'Min4x4': A minimal UPB in $\mathbb{C}^4 \otimes \mathbb{C}^4$ constructed in [3].
    • 'Pyramid': A minimal UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ constructed in [1].
    • 'QuadRes': A minimal UPB in $\mathbb{C}^d \otimes \mathbb{C}^d$ (only valid when 2d-1 is prime and d is odd) constructed in [2]. Note that you must set OPT_PAR equal to d (the local dimension) in this case.
    • 'Shifts': A minimal UPB in $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$ introduced in [1].
    • 'SixParam': The six-parameter UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ introduced in Section IV.A of [2]. Note that OPT_PAR must be a vector containing the six parameters in this case.
    • 'Tiles': A minimal UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ constructed in [1].
  • DIM: A vector containing the local dimensions of the desired UPB. In all cases, the smallest known UPB of the desired dimensionality is returned. If no unextendible product basis is known for the specified dimensions, an error is produced.
  • VERBOSE (optional, default 1): A flag (either 1 or 0) indicating whether or not a reference to a journal article that contains the UPB (or a description of how to construct the UPB) returned by this script will be displayed.

Output arguments

If only one output argument is specified (e.g., U = UPB(DIM)) then U is a matrix whose columns are the product states in the desired UPB.

If multiple output arguments are specified (e.g., [U,V,W,...] = UPB(DIM)) then the unextendible product basis is obtained by tensoring the columns of U, V, W, ... together. That is, U, V, W, ... are the local vectors in the unextendible product basis.

Examples

Generating the "Shifts" UPB

The following code returns the "Shifts" UPB [1], which is a UPB of 4 states on $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$:

>> v = UPB('Shifts')
 
v =
 
    1.0000   -0.0000   -0.0000   -0.0000
         0   -0.0000    0.0000   -0.5000
         0    0.0000   -0.5000   -0.0000
         0    0.0000    0.5000   -0.5000
         0   -0.5000   -0.0000    0.0000
         0   -0.5000    0.0000    0.5000
         0    0.5000   -0.5000    0.0000
         0    0.5000    0.5000    0.5000

Alternatively, we can request that the local vectors on each copy of $\mathbb{C}^2$ are returned, rather than the total product vectors on $\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^2$:

>> [u,v,w] = UPB('Shifts')
 
u =
 
    1.0000    0.0000    0.7071   -0.7071
         0    1.0000    0.7071    0.7071
 
 
v =
 
    1.0000   -0.7071    0.0000    0.7071
         0    0.7071    1.0000    0.7071
 
 
w =
 
    1.0000    0.7071   -0.7071    0.0000
         0    0.7071    0.7071    1.0000

Generating bound entangled states

As noted in [1], if \(\big\{|v_i\rangle\big\}\) is an unextendible product basis, then $I - \sum_i |v_i\rangle\langle v_i|$ is (up to scaling) a bound entangled state. The following code illustrates this fact in $\mathbb{C}^3 \otimes \mathbb{C}^5$ by first constructing a UPB in this space, then constructing the corresponding state, and then verifying that this state is bound entangled.

>> v = UPB([3,5]);
Generated a minimal 7-state UPB from:
N. Alon and L. Lovasz. Unextendible product bases. J. Combinatorial Theory, Ser. A, 95:169-179, 2001.
See also: http://www.njohnston.ca/2013/03/how-to-construct-minimal-upbs/
 
>> rho = eye(3*5);
>> for j = 1:7
       rho = rho - v(:,j)*v(:,j)';
   end
   rho = rho/7; % we are now done constructing the bound entangled state
 
>> IsSeparable(rho,[3,5]) % show that the state is indeed entangled
Determined to be entangled by not having a 2-copy symmetric extension. Reference:
A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri. A complete family of separability criteria.
Phys. Rev. A, 69:022308, 2004.
 
ans =
 
     0
>> IsPPT(rho,2,[3,5]) % verify that this state has positive partial transpose and is thus bound entangled
 
ans =
 
     1

Generating multipartite UPBs

Many multipartite minimal UPBs can be constructed via this script. For example, the following code generates a minimal UPB (of 6 states) in \(\mathbb{C}^2 \otimes \mathbb{C}^2 \otimes \mathbb{C}^3\):

>> [u,v,w] = UPB([2,2,3])
Generated a minimal 6-state UPB from:
K. Feng. Unextendible product bases and 1-factorization of complete graphs.
Discrete Applied Mathematics, 154:942–949, 2006.
 
u =
 
    1.0000         0    1.0000    0.7071    0.7071    0.7071
         0    1.0000         0    0.7071   -0.7071    0.7071
 
 
v =
 
    1.0000    0.7071         0         0    0.7071    1.0000
         0    0.7071    1.0000    1.0000   -0.7071         0
 
 
w =
 
    1.0000    0.7071    0.5345    0.5774         0         0
         0         0   -0.8018    0.5774    0.3162    1.0000
         0    0.7071   -0.2673   -0.5774   -0.9487         0

However, the minimum size of UPBs is still unknown in many multipartite cases – an error is returned in these cases:

>> [u,v,w,x] = UPB([2,2,3,7])
??? Error using ==> UPB at 132
No minimal UPB is currently known in the specified dimensions.

Source code

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

  1. %%  UPB  Generates an unextendible product basis
  2. %   This function may be called in several different ways:
  3. %
  4. %   U = UPB(NAME) is a matrix containing as its columns the vectors in the
  5. %   unextendible product basis specified by the string NAME. NAME must be
  6. %   one of: 'GenShifts', 'Min4x4', 'Pyramid', 'QuadRes', 'Shifts',
  7. %   'SixParam', or 'Tiles'. See the online documentation for descriptions
  8. %   of these different UPBs.
  9. %
  10. %   [U,V,W,...] = UPB(NAME) is the same as above, except in this case, the
  11. %   unextendible product basis is obtained by tensoring the columns of U,
  12. %   V, W, ... together. That is, U, V, W, ... are the local vectors in the
  13. %   unextendible product basis.
  14. %
  15. %   U = UPB(DIM) and [U,V,W,...] = UPB(DIM) are as above, except DIM is a
  16. %   vector containing the local dimensions of the requested UPB rather than
  17. %   the name of the UPB.
  18. %
  19. %   U = UPB(DIM,VERBOSE) and [U,V,W,...] = UPB(DIM,VERBOSE) are as above,
  20. %   where VERBOSE is a flag (either 1 or 0, default 1) that indicates that
  21. %   a reference for the returned UPB will or will not be displayed.
  22. %
  23. %   URL: http://www.qetlab.com/UPB
  24.  
  25. %   requires: FourierMatrix.m, IsTotallyNonsingular.m, MinUPBSize.m,
  26. %             normalize_cols.m, one_factorization.m, opt_args.m,
  27. %             opt_disp.m, perm_inv.m, PermuteSystems.m, sporth.m, Swap.m
  28. %
  29. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  30. %   package: QETLAB
  31. %   last updated: November 12, 2014
  32. %   URL: http://www.qetlab.com/UPB
  33.  
  34. function [u,varargout] = UPB(name,varargin)
  35.  
  36.     show_name = false;
  37.     given_dims = false;
  38.     revp = -1; % by default, don't permute systems around after we're done constructing the UPB
  39.     if(isnumeric(name)) % user provided dimensions, not a name, so find an appropriate UPB
  40.         % set optional argument defaults: verbose=1
  41.         [verbose] = opt_args({ 1 },varargin{:});
  42.  
  43.         given_dims = true;
  44.         np = length(name);
  45.  
  46.         % pre-load various references
  47.         refs = {'K. Feng. Unextendible product bases and 1-factorization of complete graphs. Discrete Applied Mathematics, 154:942-949, 2006.', ...
  48.                 'D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases, uncompletable product bases and bound entanglement. Commun. Math. Phys. 238, 379-410, 2003.', ...
  49.                 'C.H. Bennett, D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases and bound entanglement. Phys. Rev. Lett. 82, 5385-5388, 1999.', ...
  50.                 'T.B. Pedersen. Characteristics of unextendible product bases. Thesis, Aarhus Universitet, Datalogisk Institut, 2002.', ...
  51.                 'N. Johnston. The minimum size of qubit unextendible product bases. In Proceedings of the 8th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2013). E-print: arXiv:1302.1604 [quant-ph], 2013.', ...
  52.                 'The realm of common sense (if there is only a single party, the only UPBs are full bases of the space).', ...
  53.                 'N. Alon and L. Lovasz. Unextendible product bases. J. Combinatorial Theory, Ser. A, 95:169-179, 2001.\nSee also: http://www.njohnston.ca/2013/03/how-to-construct-minimal-upbs/', ...
  54.                 'J. Chen and N. Johnston. The minimum size of unextendible product bases in the bipartite case (and some multipartite cases). Comm. Math. Phys., 333(1):351-365, 2015.'};
  55.  
  56.         [name,revp] = sort(name);
  57.         if(np == 1)
  58.             upbp{1} = eye(name);
  59.             name = '';
  60.             ref_ind = 6;
  61.         elseif(np == 2 && min(name) <= 2)
  62.             % In this case, there are no UPBs smaller than a basis of the
  63.             % whole space.
  64.             upbp{1} = repmat(eye(name(1)),1,name(2));
  65.             upbp{2} = repmat(eye(name(2)),1,name(1));
  66.             upbp{revp(1)} = Swap(repmat(eye(name(revp(1))),1,name(revp(2))).',[1,2],[name(revp(2)),name(revp(1))],1).';
  67.  
  68.             ref_ind = 2;
  69.             name = '';
  70.         elseif(np == 2 && all(name == [3,3]))
  71.             name = 'Tiles';
  72.             ref_ind = 3;
  73.             show_name = true;
  74.         elseif(np == 2 && all(name == [4,4]))
  75.             name = 'Feng4x4';
  76.             ref_ind = 1;
  77.         elseif(np == 2 && all(name == name(1)) && mod(name(1),2) == 1 && isprime(2*name(1)-1))
  78.             varargin = {name(1)};
  79.             name = 'QuadRes';
  80.             ref_ind = 2;
  81.             show_name = true;
  82.         elseif(np == 3 && all(name == [2,2,2]))
  83.             name = 'Shifts';
  84.             ref_ind = 3;
  85.             show_name = true;
  86.         elseif(np == 3 && all(name == [2,2,3]))
  87.             name = 'Feng2x2x3';
  88.             ref_ind = 1;
  89.         elseif(np == 3 && all(name == [2,2,5]))
  90.             name = 'Feng2x2x5';
  91.             ref_ind = 1;
  92.         elseif(np == 4 && all(name == [2,2,2,2]))
  93.             name = 'Feng2x2x2x2';
  94.             ref_ind = 1;
  95.         elseif(np == 4 && all(name == [2,2,2,4]))
  96.             name = 'Feng2x2x2x4';
  97.             ref_ind = 1;
  98.         elseif(np == 5 && all(name == [2,2,2,2,5]))
  99.             name = 'Feng2x2x2x2x5';
  100.             ref_ind = 1;
  101.         elseif(np == 8 && all(name == [2,2,2,2,2,2,2,2]))
  102.             name = 'John2^8';
  103.             ref_ind = 5;
  104.         elseif(mod(np,4) == 0 && all(name == 2*ones(1,np)))
  105.             varargin = {np};
  106.             name = 'John2^4k';
  107.             ref_ind = 5;
  108.         elseif(mod(np,2) == 1 && all(name == 2*ones(1,np)))
  109.             varargin = {np};
  110.             name = 'GenShifts';
  111.             ref_ind = 2;
  112.             show_name = true;
  113.         elseif((mod(sum(name)-np,2) == 1 || sum(mod(name,2)==0) == 0) && sum(mod(name,2)==0) <= 1)
  114.             varargin = {name};
  115.             name = 'AlonLovasz';
  116.             ref_ind = 7;
  117.         elseif(name(end)-1 == sum(name(1:end-1)-1) && sum(name-1) >= 3 && mod(sum(name)-np,2) == 0)
  118.             varargin = {name};
  119.             name = 'CJBip';
  120.             ref_ind = 8;
  121.         elseif(np == 2 && all(name == [4,6]))
  122.             name = 'CJBip46';
  123.             ref_ind = 8;
  124.         elseif(np == 3 && name(1) == 2 && name(2) == 2 && mod(name(3),4) == 1)
  125.             varargin = {name(3)};
  126.             name = 'CJ4k1';
  127.             ref_ind = 8;
  128.         else
  129.             try
  130.                 min_size = MinUPBSize(name,0);
  131.             catch err
  132.                 if(strcmpi(err.identifier,'MinUPBSize:MinSizeUnknown'))
  133.                     error('UPB:MinSizeUnknown','No minimal UPB is currently known in the specified dimensions.');
  134.                 else
  135.                     rethrow(err);
  136.                 end
  137.             end
  138.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
  139.         end
  140.     end
  141.  
  142.     if(strcmpi(name,'Shifts')) % GenShifts reduces to Shifts when there are 3 parties
  143.         name = 'GenShifts';
  144.         varargin = {3};
  145.     end
  146.  
  147.     % the "Pyramid" UPB
  148.     if(strcmpi(name,'Pyramid'))
  149.         h = sqrt(1 + sqrt(5))/2;
  150.         for j = 4:-1:0 % pre-allocate
  151.             upbp{1}(:,j+1) = [cos(2*pi*j/5);sin(2*pi*j/5);h];
  152.         end
  153.         upbp{1} = 2*upbp{1}/sqrt(5+sqrt(5));
  154.         upbp{2} = upbp{1}(:,[1,3,5,2,4]);
  155.  
  156.     % the "Tiles" UPB
  157.     elseif(strcmpi(name,'Tiles'))
  158.         upbp{1}(:,5) = ones(3,1)/sqrt(3); % pre-allocate
  159.         upbp{1}(:,1) = [1;0;0];
  160.         upbp{1}(:,2) = [1;-1;0]/sqrt(2);
  161.         upbp{1}(:,3) = [0;0;1];
  162.         upbp{1}(:,4) = [0;1;-1]/sqrt(2);
  163.         upbp{2}(:,5) = ones(3,1)/sqrt(3);
  164.         upbp{2}(:,1) = [1;-1;0]/sqrt(2);
  165.         upbp{2}(:,2) = [0;0;1];
  166.         upbp{2}(:,3) = [0;1;-1]/sqrt(2);
  167.         upbp{2}(:,4) = [1;0;0];
  168.  
  169.     % the "Min4x4" UPB
  170.     elseif(strcmpi(name,'Min4x4'))
  171.         upbp{1}(:,8) = [0;0;1;0]; % pre-allocate
  172.         upbp{1}(:,1) = [1;-3;1;1]/sqrt(12);
  173.         upbp{1}(:,2) = [1;0;0;0];
  174.         upbp{1}(:,3) = [0;1;2;1]/sqrt(6);
  175.         upbp{1}(:,4) = [1;0;0;-1]/sqrt(2);
  176.         upbp{1}(:,5) = [0;1;0;0];
  177.         upbp{1}(:,6) = [3;1;-1;1]/sqrt(12);
  178.         upbp{1}(:,7) = [0;1;1;0]/sqrt(2);
  179.         upbp{2}(:,8) = [-1;1+sqrt(2);0;1]/sqrt(5+2*sqrt(2)); % pre-allocate
  180.         upbp{2}(:,1) = [0;1;-3-sqrt(2);-1-sqrt(2)]/sqrt(15+8*sqrt(2));
  181.         upbp{2}(:,2) = [1;0;0;0];
  182.         upbp{2}(:,3) = [1;0;sqrt(2)-1;1]/sqrt(5-2*sqrt(2));
  183.         upbp{2}(:,4) = [0;1;0;0];
  184.         upbp{2}(:,5) = [-1;1+sqrt(2);0;1]/sqrt(5+2*sqrt(2));
  185.         upbp{2}(:,6) = [0;0;1;0];
  186.         upbp{2}(:,7) = [1;1;1;-sqrt(2)]/sqrt(5);
  187.  
  188.     % the "QuadRes" UPB
  189.     elseif(strcmpi(name,'QuadRes'))
  190.         if(isempty(varargin))
  191.             error('UPB:InvalidArguments','When using NAME=QuadRes, you must specify a second input argument that gives the dimension of the desires UPB.')
  192.         elseif(~isprime(2*varargin{1}-1))
  193.             error('UPB:InvalidArguments','When using NAME=QuadRes, the second input argument DIM must be such that 2*DIM-1 is prime.')
  194.         elseif(mod(varargin{1},2) == 0)
  195.             error('UPB:InvalidArguments','When using NAME=QuadRes, the second input argument DIM must be odd.')
  196.         end
  197.  
  198.         p = 2*varargin{1}-1;
  199.         q = quad_residue(p);
  200.         s = setdiff(1:p-1,q);
  201.         sm = sum(exp(2i*pi*q/p));
  202.         N = max(-sm,1+sm);
  203.  
  204.         % This UPB is obtained by choosing entries from the Fourier matrix
  205.         % carefully.
  206.         F = FourierMatrix(p);
  207.         F(1,:) = sqrt(N)*F(1,:);
  208.         upbp{1} = F([1,q+1],:);
  209.         upbp{2} = F([1,mod(s(1)*q,p)+1],:);
  210.  
  211.         % normalize the output
  212.         upbp{1} = upbp{1}./repmat(sqrt(sum(abs(upbp{1}).^2,1)),varargin{1},1);
  213.         upbp{2} = upbp{2}./repmat(sqrt(sum(abs(upbp{2}).^2,1)),varargin{1},1);
  214.  
  215.     % the "SixParam" UPB (i.e., the one from Section IV.A of DMSST03)
  216.     elseif(strcmpi(name,'SixParam'))
  217.         if(isempty(varargin) || length(varargin{1}) ~= 6)
  218.             error('UPB:InvalidArguments','When using NAME=SixParam, you must specify a vector containing the six parameters [gammaA,thetaA,phiA,gammaB,thetaB,phiB].')
  219.         elseif(any(abs(sin(varargin{1}([1,2,4,5]))) < 10*eps) || any(abs(cos(varargin{1}([1,2,4,5]))) < 10*eps))
  220.             error('UPB:InvalidArguments','When using NAME=SixParam, none of the gammaA, thetaA, gammaB, and thetaB parameters can be a multiple of pi/2.')
  221.         end
  222.  
  223.         arg_cell = num2cell(varargin{1});
  224.         [gammaA,thetaA,phiA,gammaB,thetaB,phiB] = arg_cell{:}; % give more memorable names to parameters
  225.  
  226.         NA = sqrt(cos(gammaA)^2 + sin(gammaA)^2 * cos(thetaA)^2);
  227.         NB = sqrt(cos(gammaB)^2 + sin(gammaB)^2 * cos(thetaB)^2);
  228.  
  229.         upbp{1}(:,5) = [0;sin(gammaA)*cos(thetaA)*exp(1i*phiA);cos(gammaA)]/NA; % pre-allocate
  230.         upbp{1}(:,1) = [1;0;0];
  231.         upbp{1}(:,2) = [0;1;0];
  232.         upbp{1}(:,3) = [cos(thetaA);0;sin(thetaA)];
  233.         upbp{1}(:,4) = [sin(gammaA)*sin(thetaA);cos(gammaA)*exp(1i*phiA);-sin(gammaA)*cos(thetaA)];
  234.         upbp{2}(:,5) = [0;sin(gammaB)*cos(thetaB)*exp(1i*phiB);cos(gammaB)]/NB;
  235.         upbp{2}(:,1) = [0;1;0];
  236.         upbp{2}(:,2) = [sin(gammaB)*sin(thetaB);cos(gammaB)*exp(1i*phiB);-sin(gammaB)*cos(thetaB)];
  237.         upbp{2}(:,3) = [1;0;0];
  238.         upbp{2}(:,4) = [cos(thetaB);0;sin(thetaB)];
  239.  
  240.     % the "GenShifts" UPB
  241.     elseif(strcmpi(name,'GenShifts'))
  242.         if(isempty(varargin))
  243.             error('UPB:InvalidArguments','When using NAME=GenShifts, you must specify a second input argument that gives the number of parties in the desired UPB.')
  244.         elseif(mod(varargin{1},2) == 0)
  245.             error('UPB:InvalidArguments','When using NAME=GenShifts, the second input argument P must be odd.')
  246.         end
  247.  
  248.         k = (varargin{1}+1)/2;
  249.         upbp{1} = [cos((0:(1/k):(2-1/k))*pi/2);sin((0:(1/k):(2-1/k))*pi/2)];
  250.         upbp{1} = upbp{1}(:,[1,k+1:-1:2,k+2:end]);
  251.  
  252.         for j = varargin{1}-1:-1:1 % pre-allocate memory
  253.             upbp{j+1} = upbp{1}(:,[1,circshift(2:(2*k),[0,j])]);
  254.         end
  255.  
  256.     % the "Feng2x2x2x2" UPB
  257.     elseif(strcmpi(name,'Feng2x2x2x2'))
  258.         b1 = eye(2);
  259.         b2 = [1 1;1 -1]/sqrt(2);
  260.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
  261.  
  262.         upbp{1} = [b1(:,1),b1(:,2),b1(:,1),b2(:,1),b2(:,2),b2(:,1)];
  263.         upbp{2} = [b1(:,1),b2(:,1),b1(:,2),b1(:,2),b2(:,2),b1(:,1)];
  264.         upbp{3} = [b1(:,1),b2(:,1),b3(:,1),b2(:,2),b3(:,2),b1(:,2)];
  265.         upbp{4} = [b1(:,1),b2(:,1),b3(:,1),b3(:,2),b1(:,2),b2(:,2)];
  266.  
  267.     % the "John2^8" UPB
  268.     elseif(strcmpi(name,'John2^8'))
  269.         b1 = eye(2);
  270.         b2 = [1 1;1 -1]/sqrt(2);
  271.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
  272.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
  273.         b5 = [cos(pi/7) sin(pi/7);-sin(pi/7) cos(pi/7)];
  274.  
  275.         upbp{1} = [b1(:,1),b2(:,1),b2(:,2),b1(:,2),b3(:,1),b4(:,1),b4(:,1),b3(:,1),b1(:,2),b3(:,2),b4(:,2)];
  276.         upbp{2} = [b1(:,1),b1(:,2),b2(:,1),b3(:,1),b4(:,1),b4(:,1),b4(:,2),b4(:,2),b3(:,1),b2(:,2),b3(:,2)];
  277.         upbp{3} = [b1(:,1),b2(:,1),b3(:,1),b3(:,2),b1(:,2),b3(:,2),b1(:,2),b4(:,1),b3(:,1),b2(:,2),b4(:,2)];
  278.         upbp{4} = [b1(:,1),b2(:,1),b3(:,1),b3(:,1),b3(:,2),b4(:,1),b3(:,2),b2(:,2),b4(:,1),b4(:,2),b1(:,2)];
  279.         upbp{5} = [b1(:,1),b2(:,1),b1(:,2),b3(:,1),b4(:,1),b2(:,2),b5(:,1),b3(:,2),b3(:,1),b5(:,2),b4(:,2)];
  280.         upbp{6} = [b1(:,1),b2(:,1),b3(:,1),b2(:,2),b4(:,1),b4(:,2),b5(:,1),b5(:,2),b2(:,2),b1(:,2),b3(:,2)];
  281.         upbp{7} = [b1(:,1),b2(:,1),b3(:,1),b4(:,1),b2(:,2),b4(:,2),b2(:,2),b1(:,2),b3(:,2),b5(:,1),b5(:,2)];
  282.         upbp{8} = [b1(:,1),b2(:,1),b3(:,1),b4(:,1),b5(:,1),b1(:,2),b5(:,1),b3(:,2),b5(:,2),b4(:,2),b2(:,2)];
  283.  
  284.     % the "John2^4k" UPB
  285.     elseif(strcmpi(name,'John2^4k'))
  286.         if(isempty(varargin))
  287.             error('UPB:InvalidArguments','When using NAME=John2^4k, you must specify a second input argument that gives the number of parties in the desired UPB.')
  288.         elseif(mod(varargin{1},4) ~= 0 || varargin{1} <= 7)
  289.             error('UPB:InvalidArguments','When using NAME=John2^4k, the second input argument P must equal 0 (mod 4) and must be at least 8.')
  290.         end
  291.  
  292.         % Construct varargin{1}/2+2 distinct orthonormal bases of C^2.
  293.         for j = (varargin{1}/2+2):-1:1 % pre-allocate
  294.             jb = 2*pi/(2*j-1);
  295.             b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
  296.         end
  297.  
  298.         % the first 3 parties are special, so we do them separately
  299.         upbp{1} = [reshape(repmat(reshape(b(:,1,1:varargin{1}/4+1),2,varargin{1}/4+1),2,1),2,varargin{1}/2+2), reshape(repmat(reshape(b(:,2,1:varargin{1}/4+1),2,varargin{1}/4+1),2,1),2,varargin{1}/2+2)];
  300.         upbp{2} = [upbp{1}(:,1:varargin{1}/2+2),circshift(upbp{1}(:,varargin{1}/2+3:end),[0,2])];
  301.         upbp{3} = [upbp{1}(:,1:varargin{1}/2+2),circshift(upbp{1}(:,varargin{1}/2+3:end),[0,4])];
  302.  
  303.         % now do the next 2k-4 parties
  304.         if(varargin{1} > 8)
  305.             upbp{4} = [reshape(b(:,1,:),2,varargin{1}/2+2), circshift(reshape(b(:,2,:),2,varargin{1}/2+2),[0,6])];
  306.             upbp{5} = upbp{4}(:,[1:varargin{1}/2+2,reshape(fliplr(reshape(varargin{1}/2+3:varargin{1}+4,2,varargin{1}/4+1).').',1,varargin{1}/2+2)]);
  307.  
  308.             for j = 3:(varargin{1}/4-1)
  309.                 upbp{2*j} = [upbp{4}(:,1:varargin{1}/2+2),circshift(upbp{4}(:,varargin{1}/2+3:end),[0,2*(j-1)])];
  310.                 upbp{2*j+1} = [upbp{5}(:,1:varargin{1}/2+2),circshift(upbp{5}(:,varargin{1}/2+3:end),[0,2*(j-1)])];
  311.             end
  312.         end
  313.  
  314.         % Finally, do the last 2k+1 parties, which arise from finding a
  315.         % 1-factorization of the complete graph.
  316.         fac = one_factorization(varargin{1}/2+2);
  317.         resb = reshape(b,2,varargin{1}+4);
  318.         for j = 1:(varargin{1}/2+1)
  319.             upbp{varargin{1}/2+j-1} = resb(:,[fac(j,:),varargin{1}/2+2+fac(j,:)]);
  320.         end
  321.  
  322.     % the UPB from Theorem 3 of reference [8]
  323.     elseif(strcmpi(name,'CJ4k1'))
  324.         % Construct (varargin{1}+1)/2+1 distinct orthonormal bases of C^2.
  325.         for j = ((varargin{1}+1)/2+1):-1:1 % pre-allocate
  326.             jb = 2*pi/(2*j-1);
  327.             b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
  328.         end
  329.  
  330.         % the 2-dimensional parties are the same as in the Feng4m2 UPB
  331.         upbp{1} = repmat(reshape(b(:,:,1:((varargin{1}+3)/4)),2,(varargin{1}+1)/2+1),1,2);
  332.         u_ind = reshape((varargin{1}+1)/2+2:varargin{1}+3,2,(varargin{1}+3)/4);
  333.         upbp{1}(:,(varargin{1}+1)/2+2:varargin{1}+3) = upbp{1}(:,reshape(u_ind([2,1],:),1,(varargin{1}+1)/2+1));
  334.         b = reshape(b,2,varargin{1}+3);
  335.         upbp{3} = CJ_Lemma6((varargin{1}-1)/4);
  336.         upbp{2} = b(:,circshift(1:varargin{1}+3,[0,1]));
  337.         upbp{2}(:,[(varargin{1}+3)/2,varargin{1}+3]) = upbp{2}(:,[varargin{1}+3,(varargin{1}+3)/2]);        
  338.  
  339.     % the UPB from Theorem 1 of reference [8]
  340.     elseif(strcmpi(name,'CJBip'))
  341.         maxD = varargin{1}(end);
  342.         U = HollowUnitary(maxD);
  343.         upbp{np} = [eye(maxD),U];
  344.  
  345.         if(maxD >= 20 && ~IsTotallyNonsingular(U,2:maxD-2))
  346.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(maxD*2),' in this case, but their construction is complicated and not implemented by this script.']);
  347.         end
  348.  
  349.         prevdims = 0;
  350.         for j = 1:np-1
  351.             upbp{j} = CJ_Lemma5(maxD,varargin{1}(j)-1,1+prevdims);
  352.             prevdims = prevdims + varargin{1}(j)-1;
  353.         end
  354.  
  355.     % another UPB from Theorem 1 of reference [8]
  356.     elseif(strcmpi(name,'CJBip46'))
  357.         u = 1.64451358502312496885542269243;
  358.         upbp{1} = normalize_cols([3 3 3 1 1 -1 -1 2 -2 0;2 1 1 2 3 -2 0 0 -2 -1;2 1 1 1 1 5 -3 -4 2 1;2 2 3 2 2 0 2 1 3 0]);
  359.         upbp{2} = eye(6);
  360.         upbp{2} = [upbp{2}(:,1:5),normalize_cols([0 0 (u-2)/(1-u) 1 (u-1)/(u-2);u*(u-2)/(2*u-3) 0 0 (3-2*u)/(u-2) 1;1 -1-u 0 0 1/(1+u);1 1 -u 0 0;0 u-1 1 1/(1-u) 0;u 1 1 1 1])];
  361.  
  362.     % the "Feng2x2x3" UPB
  363.     elseif(strcmpi(name,'Feng2x2x3'))
  364.         b1 = eye(2);
  365.         b2 = [1 1;1 -1]/sqrt(2);
  366.  
  367.         upbp{3} = normalize_cols([1,1,2,1,0,0;0,0,-3,1,1,1;0,1,-1,-1,-3,0]);
  368.         upbp{1} = [b1(:,1),b1(:,2),b1(:,1),b2(:,1),b2(:,2),b2(:,1)];
  369.         upbp{2} = [b1(:,1),b2(:,1),b1(:,2),b1(:,2),b2(:,2),b1(:,1)];
  370.  
  371.     % the "Feng2x2x5" UPB
  372.     elseif(strcmpi(name,'Feng2x2x5'))
  373.         w = exp(pi*2i/3);
  374.         b1 = eye(2);
  375.         b2 = [1 1;1 -1]/sqrt(2);
  376.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
  377.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
  378.  
  379.         upbp{3} = normalize_cols([1,conj(w),w,0,0,0,1,0;0,w,conj(w),1,0,1,0,0;0,conj(w),0,0,0,1,w,1;0,0,w,0,1,1,conj(w),0;0,1,1,0,0,1,1,0]);
  380.         upbp{1} = [b2(:,2),b2(:,1),b1(:,2),b1(:,1),b1(:,1),b1(:,2),b2(:,1),b2(:,2)];
  381.         upbp{2} = [b4(:,2),b1(:,2),b4(:,1),b1(:,1),b3(:,1),b2(:,1),b3(:,2),b2(:,2)];
  382.  
  383.     % the "Feng4m2" UPB of Theorem 3.2 from the Feng paper
  384.     % I don't know how to compute a 1-factorization of the complement graph
  385.     % in this UPB's construction yet. The commented section below shows the
  386.     % part of the construction that I *do* know how to implement.
  387. %    elseif(strcmpi(name,'Feng4m2'))
  388. %        if(isempty(varargin))
  389. %            error('UPB:InvalidArguments','When using NAME=Feng4m2, you must specify a second input argument that gives the number of parties in the desired UPB.')
  390. %        elseif(mod(varargin{1},4) ~= 2)
  391. %            error('UPB:InvalidArguments','When using NAME=Feng4m2, the second input argument P must equal 2 (mod 4).')
  392. %        end
  393. %        
  394. %        % Construct varargin{1}/2+1 distinct orthonormal bases of C^2.
  395. %        for j = (varargin{1}/2+1):-1:1 % pre-allocate
  396. %            jb = 2*pi/(2*j-1);
  397. %            b(:,:,j) = [cos(jb) sin(jb);-sin(jb) cos(jb)];
  398. %        end
  399. %        
  400. %        upbp{1} = repmat(reshape(b(:,:,1:((varargin{1}+2)/4)),2,varargin{1}/2+1),1,2); % these are the a_j's in the proof of the theorem
  401. %        u_ind = reshape(varargin{1}/2+2:varargin{1}+2,2,(varargin{1}+2)/4);% we need to move the columns of upbp{1} around a little bit first
  402. %        upbp{1}(:,varargin{1}/2+2:varargin{1}+2) = upbp{1}(:,reshape(u_ind([2,1],:),1,varargin{1}/2+1));
  403. %        b = reshape(b,2,varargin{1}+2); % these are the vectors that will be used in all subsystems except the first
  404.  
  405.     % the "Feng4x4" UPB (Theorem 3.3(4) of Feng paper)
  406.     elseif(strcmpi(name,'Feng4x4'))
  407.         upbp{1}(:,8) = [1;-1;1;0]/sqrt(3); % pre-allocate
  408.         upbp{1}(:,1:4) = eye(4);
  409.         upbp{1}(:,5) = [0;1;1;1]/sqrt(3);
  410.         upbp{1}(:,6) = [1;0;-1;1]/sqrt(3);
  411.         upbp{1}(:,7) = [1;1;0;-1]/sqrt(3);
  412.         upbp{2}(:,8) = [0;1;1;1]/sqrt(3);
  413.         upbp{2}(:,[1,7,6,4]) = eye(4);
  414.         upbp{2}(:,5) = [1;-1;1;0]/sqrt(3);
  415.         upbp{2}(:,2) = [1;0;-1;1]/sqrt(3);
  416.         upbp{2}(:,3) = [1;1;0;-1]/sqrt(3);
  417.  
  418.     % the "Feng2x2x2x4" UPB (Theorem 3.3(3) of Feng paper)
  419.     elseif(strcmpi(name,'Feng2x2x2x4'))
  420.         % need four different bases of C^2
  421.         b1 = eye(2);
  422.         b2 = [1 1;1 -1]/sqrt(2);
  423.         b3 = [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)];
  424.         b4 = [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)];
  425.  
  426.         upbp{4}(:,8) = [1;-1;1;0]/sqrt(3); % pre-allocate
  427.         upbp{4}(:,1:4) = eye(4);
  428.         upbp{4}(:,5) = [0;1;1;1]/sqrt(3);
  429.         upbp{4}(:,6) = [1;0;-1;1]/sqrt(3);
  430.         upbp{4}(:,7) = [1;1;0;-1]/sqrt(3);
  431.  
  432.         upbp{1} = [b1 b2 b3 b4];
  433.         upbp{3} = [b1 b2 b3 b4];
  434.         upbp{2} = [b1 b2 b3 b4];
  435.         upbp{2} = upbp{2}(:,[1,3,7,5,4,2,6,8]);
  436.         upbp{3} = upbp{3}(:,[1,5,7,3,4,8,6,2]);
  437.         upbp{1} = upbp{1}(:,[1,3,7,5,8,6,2,4]);
  438.  
  439.     % the "Feng2x2x2x2x5" UPB (Theorem 3.3(5) of Feng paper)
  440.     elseif(strcmpi(name,'Feng2x2x2x2x5'))
  441.         w = exp(pi*2i/3);
  442.  
  443.         upbp{5}(:,10) = [1;conj(w);w;1;0]/2; % pre-allocate
  444.         upbp{5}(:,1:5) = eye(5);
  445.         upbp{5}(:,6) = [0;1;1;1;1]/2;
  446.         upbp{5}(:,7) = [1;0;1;w;conj(w)]/2;
  447.         upbp{5}(:,8) = [1;1;0;conj(w);w]/2;
  448.         upbp{5}(:,9) = [1;w;conj(w);0;1]/2;
  449.  
  450.         % need five different bases of C^2
  451.         upbp{1} = [eye(2), [1 1;1 -1]/sqrt(2), [cos(pi/3) sin(pi/3);-sin(pi/3) cos(pi/3)], [cos(pi/5) sin(pi/5);-sin(pi/5) cos(pi/5)], [cos(pi/7) sin(pi/7);-sin(pi/7) cos(pi/7)]];
  452.         upbp{4} = upbp{1};
  453.         upbp{3} = upbp{1};
  454.         upbp{2} = upbp{1};
  455.  
  456.         upbp{2} = upbp{2}(:,[1,3,5,7,9,10,2,4,6,8]);
  457.         upbp{3} = upbp{3}(:,[1,3,5,7,9,8,10,2,4,6]);
  458.         upbp{4} = upbp{4}(:,[1,3,5,7,9,6,8,10,2,4]);
  459.         upbp{1} = upbp{1}(:,[1,3,5,7,9,4,6,8,10,2]);
  460.     elseif(strcmpi(name,'AlonLovasz'))
  461.         dim = varargin{1};
  462.         min_size = sum(dim) - np + 1;
  463.         os = 0;
  464.  
  465.         if(mod(dim(1),2) == 0)
  466.             upbp{1} = AL_even(dim(1),min_size);
  467.         else
  468.             upbp{1} = AL_odd(dim(1),min_size,os);
  469.             os = os + (dim(1)-1)/2;
  470.         end
  471.         if(length(upbp{1}) == 1 && upbp{1} == -1)
  472.             error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
  473.         end
  474.  
  475.         for j = np:-1:2
  476.             if(mod(dim(j),2) == 0)
  477.                 upbp{j} = AL_even(dim(j),min_size);                
  478.             else
  479.                 upbp{j} = AL_odd(dim(j),min_size,os);
  480.                 os = os + (dim(j)-1)/2;
  481.             end
  482.             if(length(upbp{j}) == 1 && upbp{j} == -1)
  483.                 error('UPB:HardToConstruct',['Minimal UPBs are known to have size ',num2str(min_size),' in this case, but their construction is complicated and not implemented by this script.']);
  484.             end
  485.         end
  486.     end
  487.  
  488.     if(length(revp) == 1)
  489.         revp = 1:length(upbp);
  490.     end
  491.     revp = perm_inv(revp);
  492.     u = upbp{revp(1)};
  493.     varargout = upbp(revp(2:end));
  494.  
  495.  
  496.     % If the user just requested one output argument, tensor local vectors
  497.     % together.
  498.     if(nargout <= 1 && length(varargout) >= 1)
  499.         num_vec = size(u,2);
  500.         % Do a clever column-by-column kronecker product that avoids having to
  501.         % loop over every vector within each party. Note that this method
  502.         % relies on not using sparse matrices.
  503.         u = reshape(u,1,size(u,1),num_vec);
  504.         for k = 1:length(varargout)-1
  505.             v = reshape(varargout{k},size(varargout{k},1),1,num_vec);
  506.             u = reshape(bsxfun(@times,v,u),1,size(u,2)*size(v,1),num_vec);
  507.         end % we split the last iteration outside of the loop to reduce the number of required reshapes by 1
  508.         v = reshape(varargout{end},size(varargout{end},1),1,num_vec);
  509.         u = reshape(bsxfun(@times,u,v), size(u,2)*size(varargout{end},1),num_vec);
  510.     end
  511.  
  512.     if(given_dims)
  513.         if(show_name)
  514.             opt_disp(['Generated the ',num2str(length(upbp{1})),'-state ''',name,''' UPB from:\n',refs{ref_ind},'\n'],verbose);
  515.         else
  516.             opt_disp(['Generated a minimal ',num2str(length(upbp{1})),'-state UPB from:\n',refs{ref_ind},'\n'],verbose);
  517.         end
  518.     end
  519. end
  520.  
  521. % Computes all quadratic residues mod p (needed for the QuadRes UPB)
  522. function q = quad_residue(p)
  523.     q = unique(mod((1:floor(p/2)).^2,p));
  524. end
  525.  
  526. % Computes the vectors in an Alon-Lovasz UPB in odd dimensions.
  527. function W = AL_odd(r,c,os)
  528.  
  529.     its = 0;
  530.     lin_indep = false;
  531.  
  532.     while ~lin_indep
  533.         W = randn(r,c);
  534.  
  535.         for j = 2:c
  536.             ind = intersect(1:j-1,mod([(j-(r-1)/2-os):(j-1-os),(j+1+os):(j+(r-1)/2+os)]-1,c)+1);
  537.             tmp = null(W(:,ind)');
  538.             W(:,j) = tmp*randn(size(tmp,2),1);
  539.         end
  540.  
  541.         W = normalize_cols(W);
  542.         lin_indep = IsTotallyNonsingular(W,r);
  543.         its = its + 1;
  544.  
  545.         if(its >= 2 && ~lin_indep)
  546.             W = -1;
  547.             return
  548.         end
  549.     end
  550. end
  551.  
  552. % Computes the vectors in an Alon-Lovasz UPB in one even dimension.
  553. function W = AL_even(r,c)
  554.  
  555.     its = 0;
  556.     lin_indep = false;
  557.  
  558.     while ~lin_indep
  559.         r2 = c - r + 1;
  560.         W = randn(r,c);
  561.  
  562.         for j = 2:c
  563.             ind = setdiff(1:c,mod((j-(r2-1)/2:j+(r2-1)/2)-1,c)+1);
  564.             tmp = null(W(:,ind)');
  565.             W(:,j) = tmp*randn(size(tmp,2),1);
  566.         end
  567.  
  568.         W = normalize_cols(W);
  569.         lin_indep = IsTotallyNonsingular(W,r);
  570.         its = its + 1;
  571.  
  572.         if(its >= 2 && ~lin_indep)
  573.             W = -1;
  574.             return
  575.         end
  576.     end
  577. end
  578.  
  579. %%  CJ_LEMMA5    Construct a matrix W described by Lemma 5 of [CJ]
  580. %   This function has three required argument:
  581. %     Q: the positive integer q described by the lemma
  582. %     R: the positive integer r described by the lemma
  583. %     S: the positive integer s described by the lemma
  584. %
  585. %   W = CJ_Lemma5(Q,R,S) is a (R+1)-by-2Q matrix with columns of unit
  586. %   length that satisfy the orthogonality condition (a) and the
  587. %   nonsingularity condition (b) of Lemma 5.
  588. %
  589. %   This function has one optional argument:
  590. %     NICE (default 0): a non-negative integer that specifies how "nice"
  591. %                       the entries of W should be
  592. %
  593. %   W = CJ_Lemma5(Q,R,S,NICE) is the same as before if NICE = 0. If NICE is
  594. %   positive integer, then W will be a symbolic matrix (rather than a
  595. %   numeric matrix) whose entries are integers. In this case, the columns
  596. %   of W are no longer scaled to have unit length, but rather are scaled to
  597. %   make the entries the smallest integers possible. Lower values of NICE
  598. %   lead to smaller entries, but in general take longer to compute (and may
  599. %   not be able to be computed at all, if NICE is too small).
  600. %
  601. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
  602. %
  603. %   References:
  604. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
  605. %        Bases in the Bipartite Case (and Some Multipartite Cases).
  606. %        Preprint, 2012.
  607.  
  608. %   requires: IsTotallyNonsingular.m, normalize_cols.m
  609. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  610. %   version: 1.00
  611. %   last updated: December 19, 2012
  612.  
  613. function W = CJ_Lemma5(q,r,s,varargin)
  614.  
  615.     if(q < r + s)
  616.         error('CJ_Lemma5:InvalidArgs','The arguments must satisfy q >= r + s.');
  617.     end
  618.  
  619.     if(nargin == 3)
  620.         nice = 0;
  621.     else
  622.         nice = varargin{1};
  623.     end
  624.  
  625.     its = 0;
  626.     satisfy_cond_b = false;
  627.     while ~satisfy_cond_b
  628.         % Construct a matrix W that satisfies condition (a) of Lemma 4.
  629.         if(nice == 0)
  630.             W = randn(r+1,2*q);
  631.             W(:,q+1:(2*q)) = normalize_cols(W(:,q+1:(2*q)));
  632.         else
  633.             W = sym((1-2*floor(2*rand(r+1,2*q))).*(1+floor(nice*rand(r+1,2*q))));
  634.         end
  635.         lin_depen = 0;
  636.  
  637.         % Now generate the remaining conditions according to condition (b).
  638.         for j = 0:q-1
  639.             tmp_null = null(W(:,q+1+mod(j+(s:s+r-1),q)).');
  640.             if(size(tmp_null,2) > 1)
  641.                 % Uh-oh, two of the randomly-generated columns are linearly
  642.                 % dependent: try again!
  643.                 lin_depen = 1;
  644.                 break;
  645.             end
  646.             W(:,j+1) = tmp_null;        
  647.  
  648.             % Get rid of the denominators in the new column, if the user wanted
  649.             % all integers.
  650.             if(nice > 0) 
  651.                 [~,cd] = numden(W(1,j+1));
  652.                 for k=2:size(W,1)
  653.                     [~,den] = numden(W(k,j+1));
  654.                     cd = lcm(cd,den);
  655.                 end
  656.  
  657.                 W(:,j+1) = cd*W(:,j+1);
  658.             end
  659.         end
  660.  
  661.         % With probability 1, the matrix W also satisfies condition (b) of
  662.         % Lemma 4. Nevertheless, we should make sure that this is the case, and
  663.         % if it isn't, try again!
  664.         if(lin_depen == 0)
  665.             W = fliplr(W);
  666.             satisfy_cond_b = IsTotallyNonsingular(double(W),r+1);
  667.         end
  668.  
  669.         its = its + 1;
  670.         if(its == 10 && nice ~= 0 && ~satisfy_cond_b)
  671.             warning('CJ_Lemma5:LinearDependence', 'Tried (and failed) 10 times to find a matrix satisfying all imposed conditions. You may have better luck if you increase (or omit) the NICE argument.');
  672.         end
  673.     end
  674. end
  675.  
  676. %%  CJ_LEMMA6    Construct a matrix W described by Lemma 6 of [CJ]
  677. %   This function has one required argument:
  678. %     K: the positive integer k described by the lemma
  679. %
  680. %   W = CJ_Lemma6(K) is a (4K+1)-by-(4K+4) matrix with columns of unit
  681. %   length that satisfy the orthogonality conditions (i), (ii) and (iii)
  682. %   and the nonsingularity condition (iv) of Lemma 6.
  683. %
  684. %   This function has one optional argument:
  685. %     NICE (default 0): a non-negative integer that specifies how "nice"
  686. %                       the entries of W should be
  687. %
  688. %   W = CJ_Lemma6(K,NICE) is the same as before if NICE = 0. If NICE is a
  689. %   positive integer, then W will be a symbolic matrix (rather than a
  690. %   numeric matrix) whose entries are integers. In this case, the columns
  691. %   of W are no longer scaled to have unit length, but rather are scaled to
  692. %   make the entries the smallest integers possible. Lower values of NICE
  693. %   lead to smaller entries, but in general take longer to compute (and may
  694. %   not be able to be computed at all, if NICE is too small).
  695. %
  696. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
  697. %
  698. %   References:
  699. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
  700. %        Bases in the Bipartite Case (and Some Multipartite Cases).
  701. %        Preprint, 2012.
  702.  
  703. %   requires: IsTotallyNonsingular.m, normalize_cols.m
  704. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  705. %   version: 1.00
  706. %   last updated: December 21, 2012
  707.  
  708. function W = CJ_Lemma6(k,varargin)
  709.  
  710.     d = 4*k + 1;
  711.  
  712.     if(nargin == 1)
  713.         nice = 0;
  714.     else
  715.         nice = varargin{1};
  716.     end
  717.  
  718.     its = 0;
  719.     satisfy_cond_b = false;
  720.     while ~satisfy_cond_b
  721.         % Construct a matrix W that satisfies conditions (i), (ii) and (iii) of
  722.         % Lemma 6.
  723.         if(nice == 0)
  724.             W = randn(d,d+3);
  725.             W(:,1:2) = normalize_cols(W(:,1:2));
  726.         else
  727.             W = sym((1-2*floor(2*rand(d,d+3))).*(1+floor(nice*rand(d,d+3))));
  728.         end
  729.  
  730.         for j = 2:d+2
  731.             hlf = (d+3)/2;
  732.             if(j < hlf)
  733.                 ind = setdiff(0:j-2,mod(j+1,hlf))+1;
  734.             else
  735.                 ind = setdiff(0:j-1,[j-hlf,hlf+mod(j-1,hlf),hlf+mod(j+1,hlf)])+1;
  736.             end
  737.             tmp = null(W(:,ind).');
  738.             W(:,j+1) = tmp(:,1);
  739.  
  740.             % Get rid of the denominators in the new column, if the user wanted
  741.             % all integers.
  742.             if(nice > 0) 
  743.                 [~,cd] = numden(W(1,j+1));
  744.                 for k=2:size(W,1)
  745.                     [~,den] = numden(W(k,j+1));
  746.                     cd = lcm(cd,den);
  747.                 end
  748.  
  749.                 W(:,j+1) = cd*W(:,j+1);
  750.             end
  751.         end
  752.  
  753.         % With probability 1, the matrix W also satisfies condition (iv) of
  754.         % Lemma 6. Nevertheless, we should make sure that this is the case, and
  755.         % if it isn't, try again!
  756.         satisfy_cond_b = IsTotallyNonsingular(double(W),d);
  757.  
  758.         its = its + 1;
  759.         if(its == 10 && nice ~= 0 && ~satisfy_cond_b)
  760.             warning('CJ_Lemma6:LinearDependence', 'Tried (and failed) 10 times to find a matrix satisfying all imposed conditions. You may have better luck if you increase (or omit) the NICE argument.');
  761.         end
  762.     end
  763. end
  764.  
  765. %%  HOLLOWUNITARY    Produces a hollow unitary matrix
  766. %   This function has one required argument:
  767. %     DIM: a positive integer not equal to 1 or 3
  768. %
  769. %   U = HollowUnitary(DIM) is a unitary matrix with all of its diagonal
  770. %   entries equal to 0 and all of its off-diagonal entries non-zero. Such a
  771. %   matrix exists if and only if DIM = 2 or DIM >= 4.
  772. %
  773. %   URL: http://www.njohnston.ca/publications/minimum-upbs/code/
  774. %
  775. %   References:
  776. %   [CJ] J. Chen and N. Johnston. The Minimum Size of Unextendible Product
  777. %        Bases in the Bipartite Case (and Some Multipartite Cases).
  778. %        Preprint, 2012.
  779.  
  780. %   requires: FourierMatrix.m
  781. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  782. %   version: 1.00
  783. %   last updated: December 19, 2012
  784.  
  785. function U = HollowUnitary(dim)
  786.  
  787.     if(dim == 3 || dim <= 1)
  788.         error('HollowUnitary:InvalidDim','Hollow unitaries only exist when DIM = 2 or DIM >= 4.');
  789.     end
  790.  
  791.     w2 = exp(2i*pi/(dim-2));
  792.     F = FourierMatrix(dim-1);
  793.  
  794.     U2 = zeros(dim-1);
  795.     for j = 1:dim-2
  796.         U2 = U2 + w2^j*F(:,j+1)*F(:,j+1)';
  797.     end
  798.  
  799.     U = [0,ones(1,dim-1)/sqrt(dim-1);ones(dim-1,1)/sqrt(dim-1),U2];
  800. end

References

  1. 1.0 1.1 1.2 1.3 1.4 1.5 C.H. Bennett, D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases and bound entanglement. Phys. Rev. Lett. 82, 5385–5388, 1999. E-print: arXiv:quant-ph/9808030
  2. 2.0 2.1 2.2 D.P. DiVincenzo, T. Mor, P.W. Shor, J.A. Smolin, and B.M. Terhal. Unextendible product bases, uncompletable product bases and bound entanglement. Commun. Math. Phys. 238, 379–410, 2003. E-print: arXiv:quant-ph/9908070
  3. T.B. Pedersen. Characteristics of unextendible product bases. Master's Thesis, Aarhus Universitet, Datalogisk Institut, 2002.