Logo Search packages:      
Sourcecode: semidef-oct version File versions  Download package

mtrx_nrm.m

function [x,Q] = mtrx_nrm(A,q)
%
% [x,Q] = mtrx_nrm(A,q);
%
% solves the matrix norm minimization problem 
%
% (primal) minimize \| A_0 + x_1*A_1 + ... + x_k*A_k \|
%
% (dual)   maximize    - Tr A_0*Q
%          subject to  Tr A_i*Q = 0
%                      sum of singular values of Q is less than one.
% 
% The matrices A_i are of size p x q.
%
% The problem is formulated as a positive definite program:
% 
%    minimize    t 
% 
%                [ t*I     A(x) ] 
%    subject to  [              ]  >= 0 
%                [ A(x)^T  t*I  ]  
%
% with A(x) = A_0 + x_1*A_1 + ... + x_k*A_k.
%  
% Input arguments:
% A:        p times (k+1)*q;  [A_0 A_1 ... A_k].
% q:        number of columns of A_i.
%
% Output arguments:
% x:        k-vector. primal solution. 
% Q:        (p times q)-matrix. dual solution.
%
% The relative duality gap is less than 0.1 percent.

% dimensions 
p = size(A,1);  k = size(A,2)/q - 1;

#printf("mtrx_nrm: p= %d, q=%d, k=%d\n",p,q,k);

if (k < 1) | (k - fix(k) > 1e-10),   
   error('Number of columns of A must be a multiple of q.');
end;
n = p+q;  m = k+1; 

% F_i = ,  i=1,...,k+1;  F_{k+2} = I
A = reshape(A,p*q,k+1);
F = zeros(n,n*(m+1));
for i=1:k+1
   F(:,(i-1)*n+[1:n]) = [ zeros(p,p)            reshape(A(:,i),p,q)
                          reshape(A(:,i),p,q)'  zeros(q,q)          ];
end;
F(:,(k+1)*n+[1:n]) = eye(n);
F = reshape(F,n*n,m+1);
c = ;
blck_szs = n;

#printf("matrix_nrm: F=\n");
#F

% primal initial solution: 
%  x0 minimizes \| A(x) \|_F;  t = 0.001 + norm of A(x0)
x0 = -A(:,2:k+1) \ A(:,1);
x0 = ;

#printf("mtrx_nrm: x0=\n");
#x0

% dual initial solution:
% Z_{12} = U*S*V' is residual of ls solution, scaled to have
%    nuclear norm less than 1/2
% Z_{11} = U*S*U' + gamma*I,
% Z_{22} = V*S*V' + gamma*I
% with gamma = (1 - 2*nucnorm)/n  (makes trace Z0 equal to 1)

Q = -reshape(A*[1; x0(1:k)],p,q);  % residual of ls solution
[U,S,V] = svd(Q);                   
S = diag(S(1:min(p,q),1:min(p,q))); % diagonal of S
U = U(:,1:min(p,q));
V = V(:,1:min(p,q));
nucnorm = sum(S);           
if nucnorm > 0.499,   % scale to make nuclear norm less than 0.499
   Q = (0.499/nucnorm)*Q;  
   S = (0.499/nucnorm)*S;  
   nucnorm = 0.499;
end;
Z0 = [ (U .* S(:,ones(1,p))') * U'    Q 
        Q'                           (V .* S(:,ones(1,q))') * V' ] ...
     + ((1 - 2*nucnorm)/n) * eye(n)  ;
Z0 = Z0(:);

nu = 20.0;  abstol = 1e-10;  reltol = 1e-3;  tv = 0.0;  maxiters = 100;
[x,Z,ul,info,time] = sp(F,blck_szs,c,x0,Z0,nu,abstol,reltol,tv, ...  
                        maxiters);
if time(3) == maxiters,
   error('Maximum number of iterations exceeded.');
end;

% Q = (1,2)-block of Z
Z = reshape(Z,n,n); 
Q = 2*Z(1:p,p+[1:q]);

Generated by  Doxygen 1.6.0   Back to index