UPB

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', 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 . 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 .
  • 'Pyramid': A minimal UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ constructed in .
  • '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 . 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 .
  • 'SixParam': The six-parameter UPB in $\mathbb{C}^3 \otimes \mathbb{C}^3$ introduced in Section IV.A of . 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 .
• 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 , 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 , 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

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

34. function [u,varargout] = UPB(name,varargin)

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{:});

43.         given_dims = true;

44.         np = length(name);

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.'};

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).';

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

142.     if(strcmpi(name,'Shifts')) % GenShifts reduces to Shifts when there are 3 parties

143.         name = 'GenShifts';

144.         varargin = {3};

145.     end

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]);

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];

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);

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

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);

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],:);

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);

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

223.         arg_cell = num2cell(varargin{1});

224.         [gammaA,thetaA,phiA,gammaB,thetaB,phiB] = arg_cell{:}; % give more memorable names to parameters

226.         NA = sqrt(cos(gammaA)^2 + sin(gammaA)^2 * cos(thetaA)^2);

227.         NB = sqrt(cos(gammaB)^2 + sin(gammaB)^2 * cos(thetaB)^2);

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)];

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

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]);

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

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)];

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)];

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)];

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)];

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

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

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])];

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)]);

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

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

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

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]);        

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];

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

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

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])];

362.     % the "Feng2x2x3" UPB

363.     elseif(strcmpi(name,'Feng2x2x3'))

364.         b1 = eye(2);

365.         b2 = [1 1;1 -1]/sqrt(2);

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)];

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)];

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)];

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

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);

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)];

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);

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]);

439.     % the "Feng2x2x2x2x5" UPB (Theorem 3.3(5) of Feng paper)

440.     elseif(strcmpi(name,'Feng2x2x2x2x5'))

441.         w = exp(pi*2i/3);

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;

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

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;

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

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

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));

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

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

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

526. % Computes the vectors in an Alon-Lovasz UPB in odd dimensions.

527. function W = AL_odd(r,c,os)

529.     its = 0;

530.     lin_indep = false;

532.     while ~lin_indep

533.         W = randn(r,c);

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

541.         W = normalize_cols(W);

542.         lin_indep = IsTotallyNonsingular(W,r);

543.         its = its + 1;

545.         if(its >= 2 && ~lin_indep)

546.             W = -1;

547.             return

548.         end

549.     end

550. end

552. % Computes the vectors in an Alon-Lovasz UPB in one even dimension.

553. function W = AL_even(r,c)

555.     its = 0;

556.     lin_indep = false;

558.     while ~lin_indep

559.         r2 = c - r + 1;

560.         W = randn(r,c);

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

568.         W = normalize_cols(W);

569.         lin_indep = IsTotallyNonsingular(W,r);

570.         its = its + 1;

572.         if(its >= 2 && ~lin_indep)

573.             W = -1;

574.             return

575.         end

576.     end

577. end

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.

608. %   requires: IsTotallyNonsingular.m, normalize_cols.m

609. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)

610. %   version: 1.00

611. %   last updated: December 19, 2012

613. function W = CJ_Lemma5(q,r,s,varargin)

615.     if(q < r + s)

616.         error('CJ_Lemma5:InvalidArgs','The arguments must satisfy q >= r + s.');

617.     end

619.     if(nargin == 3)

620.         nice = 0;

621.     else

622.         nice = varargin{1};

623.     end

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;

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;        

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

657.                 W(:,j+1) = cd*W(:,j+1);

658.             end

659.         end

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

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

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.

703. %   requires: IsTotallyNonsingular.m, normalize_cols.m

704. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)

705. %   version: 1.00

706. %   last updated: December 21, 2012

708. function W = CJ_Lemma6(k,varargin)

710.     d = 4*k + 1;

712.     if(nargin == 1)

713.         nice = 0;

714.     else

715.         nice = varargin{1};

716.     end

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

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);

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

749.                 W(:,j+1) = cd*W(:,j+1);

750.             end

751.         end

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);

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

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.

780. %   requires: FourierMatrix.m

781. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)

782. %   version: 1.00

783. %   last updated: December 19, 2012

785. function U = HollowUnitary(dim)

787.     if(dim == 3 || dim <= 1)

788.         error('HollowUnitary:InvalidDim','Hollow unitaries only exist when DIM = 2 or DIM >= 4.');

789.     end

791.     w2 = exp(2i*pi/(dim-2));

792.     F = FourierMatrix(dim-1);

794.     U2 = zeros(dim-1);

795.     for j = 1:dim-2

796.         U2 = U2 + w2^j*F(:,j+1)*F(:,j+1)';

797.     end

799.     U = [0,ones(1,dim-1)/sqrt(dim-1);ones(dim-1,1)/sqrt(dim-1),U2];

800. end