In this paper, we use an algorithm due to G.N. Watson for finding a lower
bound for the norm of a Hadamard multiplier to find a Haagerup factorization
which leads to an upper bound and combine the two into an algorithm for
finding the norm of a Hadamard multiplier. Although we have not proved that
the algorithm always converges, it appears to always work and is, in any
case, self checking in that the program computes upper and lower bounds: if
these are (nearly) equal, the norm has been computed, if they are not, it
has not.
A Matlab implementation of the algorithm is appended below. It is likely
that this will be most helpful if you are reading the paper to get the
notation and strategy.
Download the paper in PDF (Acrobat) format.
Download the paper in ps (PostScript) format.
Download the paper in DVI format.
If you are having trouble downloading, send email to cowen@math.purdue.edu
function [kbl,kbu,S,R,V,W,x,y]=FindKB(B,Kstop,U0)
% [kbl,kbu,S,R,V,W,x,y]=FindKB(B,Kstop,U0) finds upper and lower bounds for
% the Hadamard multiplier norm of the real square matrix B. The
% lower estimate, kbl, is found using Watson's algorithm with
% Kstop iterations and the initial (unitary) matrix U0. If the
% optional argument U0 is not present, a random unitary
% is generated for the starting matrix. The upper bound, kbu,
% is found from a Haagerup factorization of B=S'*R computed using the
% data from the Watson algorithm. The upper bound algorithm
% will probably fail if there is a proper submatrix of B with
% the same Hadamard multiplier norm as B. V*W' is the computed
% value for the optimal unitary, i.e. a unitary satisfying
% K_B = norm(B.*(V*W')). The optimal rank one matrix is x*y', that is,
% it is a rank one satisfying trace-norm(B.*(x*y'))=K_B.
%
% declaring
% global N E M G D C cR cS
% will enable one to check final values of all process variables
%
%
% Set the starting unitary
%
W=eye(size(B));
if nargin==2,[V,a]=qr(rand(size(B)));else V=U0;end;
%
% Watson's algorithm: For Kstop iterations, successively,
% for U=V*W', find rank one G, trace(G' * (B.*U) ) = sp norm(B.*U) and
% for G=N(:,1)*M(:,1)', find unitary U, trace(U' * (B.*G)) = tr norm(B.*G)
%
for k=1:Kstop
[N E M]=svd(B.*(V*W'));
[V D W]=svd(B.*(N(:,1)*M(:,1)'));
end;
%
% Use convergents to find lower and upper bounds for K_B
%
[N E M]=svd(B.*(V*W'));
kbl=E(1,1);
x=N(:,1);
y=M(:,1);
%
% Check X, Y for invertibility: necessary for factoriztion
%
if min(abs(x))<10*eps | min(abs(y))<10*eps,
kbu=['infinity'];R=[];S=[];
disp('X or Y not invertible, no factorization found!')
return;
end;
%
% find factorization of B: B= S' * R, then
% find upper estimate from column norms of R and S
% Note that because S'*R=B, this is an upper bound
% no matter what happened in computing kbl
%
X=diag(x);
Y=diag(y);
C=diag(sqrt(diag(D)));
S=C*V'*inv(X);
R=inv(S')*B;
cS=zeros(size(x));
cR=cS;
for j=1:length(x)
cR(j)=norm(R(:,j));
cS(j)=norm(S(:,j));
end;
kbu=max(cR)*max(cS);
The following is a Matlab algorithm based on Theorem 1 that builds a
factorization and upper bound estimate in the Hadamard reducible case.
In a typical
case, one would apply FindKB to a matrix. If the upper and lower
estimates of the Hadamard multiplier norm are significantly different, then
the non-zero components of x and y indicate the rows and columns,
respectively, of the matrix that comprise the important submatrix.
The original matrix must be permuted so as to put that submatrix in the
upper left corner. For example, the Matlab command
function [newkbu,SS,RR,cSS,cRR,Z]=Multistep(BB,size11,Kstop,Tol)
% [newkbu,SS,RR,cSS,cRR]=Multistep(BB,size11,Kstop,Tol) finds an upper
% bound for the Hadamard multiplier norm of the real square matrix BB
% under the assumption that the multiplier norm of the upper left
% size11 by size11 corner of BB is the same as the multiplier norm
% of BB using the strategy of Theorem 1. Kstop is the number of
% iterations used in the computation of multiplier norms of the
% submatrices B11 and Z by FindKB. If the optional argument Tol
% is not present, Tol is set to 10^{-9}. If the upper and lower bounds
% of the Hadamard multiplier norm of the upper corner do not agree within
% Tol, an error message is generated and the program halted. If the
% condition of Theorem 1 for the multiplier norm of the submatrix being
% the same as the multiplier norm of BB, a warning message is generated
% and the program halted. Output is the multiplier norm of BB, the
% computed optimal factorization, the column norms of the optimal
% factorization (as a check against the theory), and Z. Z is useful if
% it is needed as input in another Multistep stage.
%
if nargin==3,Tol=1e-9; end;
nn=max(size(BB));
kk=nn-size11;
%
% Partition matrix
%
B11=BB(1:size11,1:size11);
B12=BB(1:size11,(size11+1):nn);
B21=BB((size11+1):nn,1:size11);
B22=BB((size11+1):nn,(size11+1):nn);
%
% factor corner and check accuracy of multiplier norm estimate for corner
%
[kbl,kbu,S11,R11,V,W,x,y]=FindKB(B11,Kstop);
if abs(kbu-kbl)>Tol,
disp('K_B11 not found accurately! Increase iterations or tolerance?')
return;
end;
%
% Compute upper parts of S and R and the matrix Z
%
S12=inv(R11')*B21';
R12=inv(S11')*B12;
PHI=zeros(kk,kk);
PSI=PHI;
for j=1:kk
if kbl-norm(S12(:,j))^2<0,
disp(['column ', int2str(j), ' of S12 is too big! Wrong submatrix?'])
return;
end;
if kbl-norm(R12(:,j))^2<0,
disp(['column ', int2str(j), ' of R12 is too big! Wrong submatrix?'])
return;
end;
PHI(j,j)=1/sqrt(kbl-norm(S12(:,j))^2);
PSI(j,j)=1/sqrt(kbl-norm(R12(:,j))^2);
end;
Z=PHI*(B22-B21*inv(B11)*B12)*PSI;
%
% Find the multiplier norm of Z and check condition of Theorem 1
%
[kbzl,kbzu,Sz,Rz,Vz,Wz,xz,yz]=FindKB(Z,Kstop);
if kbzu=='infinity',
disp('Upper bound for K_Z not computed! Check for subsubmatrix!')
disp(sprintf('a lower bound for K_Z is %6.4f\n', kbzl))
return;
end;
if kbzu>1,
disp('K_Z >1! K_B not equal to K_B11!')
return;
end;
%
% Compute final factorization of BB and multiplier norm of BB
%
S22=Sz*inv(PHI);
R22=Rz*inv(PSI);
SS=[S11 S12;zeros(kk,size11) S22];
RR=[R11 R12;zeros(kk,size11) R22];
cSS=zeros(1,nn);
cRR=cSS;
for j=1:nn
cRR(j)=norm(RR(:,j));
cSS(j)=norm(SS(:,j));
end;
newkbu=max(cRR)*max(cSS);