function [X,y] = CorNewton3Mex_Wnorm(G,b,I,J,tau,W) %%%%%%%%%% This code is designed to solve %%%%%% %% min 0.5||W^1/2*(X - G)*W^1/2 ||^2 %% s.t. X_ij = b_ij for (i,j) in (I,J) %% X >= tau*I X is symmetric and positive semi-definite (tau can be zero) %% %%%%%%%%%%%%%%%% Based on the algorithm in %%%%%%% %%% ``A Quadratically Convergent Newton Method for %%% Computing the Nearest Correlation Matrix'' %%% SIAM J. Matrix Anal. Appl. 28 (2006) 360--385. %%%%%%%%%%%%% By Houduo Qi and Defeng Sun %%%%%%%%%% %%%%%%%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Input % G the given symmetric correlation matrix % b the right hand side of equality constraints % I row indices of the fixed elements % J column indices of the fixed elements % tau the lower bound of the smallest eigenvalue of X* % W the weight matrix for G % % Output % X the optimal primal solution % y the optimal dual solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Last modified by Yan Gao and Defeng Sun on September 8, 2009 t0 = clock; n = length(G); k = length(b); c = ones(k,1); fprintf('\n ******************************************************** \n') fprintf( ' Semismooth Newton-CG Method for the W-norm case ') fprintf('\n ******************************************************** \n') fprintf('\n The information of this problem is as follows: \n') fprintf(' Dim. of sdp constr = %d \n',n) fprintf(' Num. of equality constr = %d \n',k) % set parameters Iter_Whole = 200; Iter_inner = 20; % maximum number of Line Search in Newton method maxit = 200; % maximum number of iterations in PCG error_tol = 1.0e-6; % termination tolerance tol = 1.0e-2; % relative accuracy for CG sigma = 1.0e-4; % tolerance in the line search of the Newton method eig_time = 0; const_sparse = 2; % check if sparse form for X should be explioted if nargin==5 W = eye(n); elseif nargin==4 tau = 0; W = eye(n); end G = (G + G')/2; W = real(W); W = (W + W')/2; w = ones(n,1); W_inv = eye(n); % reset G and b0 G0 = G; G = G - tau*eye(n); Ind = find(I==J); b0 = b; b0(Ind) = b0(Ind) - tau; % compute W^(1/2) & W^(-1/2) %% test if W is a diagonal matrix diag_if = norm(W-diag(diag(W)),'fro'); if diag_if <1.0e-12 diag_index =1; w = diag(W); if min(w) <= 0 fprintf('Warning: W is not positive definite!') return; end else diag_index =0; end t1 = clock; if diag_index == 0 % not a diagonal weight matrix [P,D] = eig(W); eig_time = eig_time + etime(clock,t1); lambda = diag(D); P = real(P); lambda = real(lambda); if lambda(1)error_tol && k1 1.0e-8 ) y = x0 + 0.5^k_inner*d; X = zeros(n,n); for i = 1:k X(I(i),J(i)) = y(i); end X = 0.5*(X + X'); if diag_index ==0 if k <= const_sparse *n X = W_half_inv*sparse(X); X = X * W_half_inv; else X = W_half_inv*X*W_half_inv; end else for i=1:n X(i,:) = X(i,:)/w_half(i); end for j=1:n X(:,j) = X(:,j)/w_half(j); end end X = G + X; X = (X + X')/2; t1 = clock; [P,D] = eig(X); eig_time = eig_time + etime(clock,t1); lambda = diag(D); P = real(P); lambda = real(lambda); if lambda(1)= num-1 && f_hist(1)-f_hist(num) < 1.0e-7 ) fprintf('\n Progress is too slow! :( ') break; end end % end for outer loop % Optimal solution X* Ip = find(lambda>0); r = length(Ip); if (r==0) X = zeros(n,n); elseif (r==n) X = X; elseif (r<=n/2) lambda1 = lambda(Ip); lambda1 = lambda1.^0.5; P1 = P(:,Ip); if r >1 P1 = P1*sparse(diag(lambda1)); X = P1*P1'; % Optimal solution X* else X = lambda1^2*P1*P1'; end else lambda2 = -lambda(r+1:n); lambda2 = lambda2.^0.5; P2 = P(:,r+1:n); P2 = P2*sparse(diag(lambda2)); X = X + P2*P2'; % Optimal solution X* end X = (X+X')/2; % rank of X* and corresponding multiplier Z* r_X = r; r_Z = length(find(abs(max(0,lambda)-lambda) > 1.0e-18)); % optimal primal and dual objective values dual_val = val_G - f; prim_val = sum(sum((X-G).*(X-G)))/2; % recover X* to the original one if diag_index ==0 if k <= const_sparse *n X = W_half_inv*sparse(X); X = X * W_half_inv; else X = W_half_inv*X*W_half_inv; end ; else for i=1:n X(i,:) = X(i,:)/w_half(i); end for j=1:n X(:,j) = X(:,j)/w_half(j); end end X = X + tau*eye(n); X = (X + X')/2; G0 = X - G0; G0 = (G0 + G0')/2; prim_val0 = sum((sum(G0.*G0)))/2; time_used = etime(clock,t0); fprintf('\n') fprintf('\n ================ Final Information ================= \n'); fprintf(' Total number of iterations = %2.0f \n',k1); fprintf(' Number of func. evaluations = %2.0f \n',f_eval) fprintf(' Number of CG Iterations = %2.0f \n',iter_cg) fprintf(' Primal objective value = %d \n', prim_val0) fprintf('\n The following results are for the transformed problem \n') fprintf(' Primal objective value = %d \n', prim_val) fprintf(' Dual objective value = %d \n', dual_val) fprintf(' Norm of gradient = %3.2e \n', norm_b) fprintf(' Rank of X* - (tau * I) === %2.0d \n', r_X) fprintf(' Rank of optimal multiplier Z* === %2.0d \n', r_Z) fprintf(' Computing time for preconditioner = %3.1f \n',prec_time) fprintf(' Computing time for CG Iterations = %3.1f \n',pcg_time) fprintf(' Computing time for eigen-decom = %3.1f \n',eig_time) fprintf(' Total Computing time (secs) = %3.1f \n',time_used) fprintf(' ====================================================== \n'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% end of the main program %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% ************************************** %% ******** All Sub-routines *********** %% ************************************** %%% To change the format of time function [h,m,s] = time(t) t = round(t); h = floor(t/3600); m = floor(rem(t,3600)/60); s = rem(rem(t,60),60); return %%% End of time.m %%% To generate the essential part of the first-order difference of d function Omega12 = omega_mat(lambda) % We compute omega only for 1<=|idx|<=n-1 n = length(lambda); idx.idp = find(lambda>0); idx.idm = setdiff([1:n],idx.idp); r = length(idx.idp); if ~isempty(idx.idp) if (r == n) Omega12 = ones(n,n); else s = n-r; dp = lambda(1:r); dn = lambda(r+1:n); Omega12 = (dp*ones(1,s))./(abs(dp)*ones(1,s) + ones(r,1)*abs(dn')); %Omega12 = max(1e-15,Omega12); %Omega = [ones(r) Omega12; Omega12' zeros(s)]; end else Omega12 = []; end return %%% End of omega_mat.m %%% To generate F(y) function [f,Fy] = gradient(y,I,J,lambda,WP,b0) n = length(WP); k = length(y); f = 0.0; Fy = zeros(k,1); const_sparse = 2; I1 = find(lambda>0); r = length(I1); if (r>0) % if (r1 WP1 = WP1*sparse(diag(lambda1)); else WP1 = lambda1*WP1; end WP1T = WP1'; if (k<=const_sparse*n) %% sparse form i=1; while (i<=k) Fy(i) = WP1(I(i),:)*WP1T(:,J(i)); i=i+1; end else %% dense form WP = WP1*WP1'; i=1; while (i<=k) Fy(i) = WP(I(i),J(i)); i=i+1; end end % else % n/2<=r<=n % lambda2 = -lambda(r+1:n); % f = lambda'*lambda - lambda2'*lambda2; % % lambda2 = lambda2.^0.5; % WP2 = WP(:, r+1:n); % WP2 = WP2*sparse(diag(lambda2)); % WP2T = WP2'; % % if (k=n/2, use a complementary formula. %H = ((E-Omega).*(WP'*sparse(Z)*WP))*WP'; H2 = P2'*sparse(Z); Omega12 = ones(r,s)- Omega12; Omega12 = Omega12.*((H2*P1)'); H = [Omega12*P2'; Omega12'*P1' + (H2*P2)*P2']; if diag_index == 1 for i=1:n Z(i,:) = Z(i,:)/w(i); end i=1; while (i<=k) Ax(i) = Z(I(i),J(i))/w(J(i)) - WP(I(i),:)*H(:,J(i)); Ax(i) = Ax(i) + 1.0e-10*x(i); i=i+1; end else Z = W_inv * sparse(Z); i=1; while (i<=k) Ax(i) = Z(I(i),:)*W_inv(:,J(i)) - WP(I(i),:)*H(:,J(i)); Ax(i) = Ax(i) + 1.0e-10*x(i); i=i+1; end end end else %dense form %Z = full(Z); to use the full form % dense form if (r<3*n/4) %H = WP*(Omega.*(WP'*Z*WP))*WP'; H1 = P1'*Z; Omega12 = Omega12.*(H1*P2); H = P1*((H1*P1)*P1'+ 2.0*Omega12*P2'); H = (H + H')/2; i = 1; while (i<=k) Ax(i) = H(I(i),J(i)); Ax(i) = Ax(i) + 1.0e-10*x(i); i = i+1; end else % if r>=3*n/4, use a complementary formula. %H = -WP*( (E-Omega).*(WP'*Z*P) )*WP'; H2 = P2'*Z; Omega12 = ones(r,s) - Omega12; Omega12 = Omega12.*(H2*P1)'; H = P2*( 2.0*(Omega12'*P1') + (H2*P2)*P2'); H = (H + H')/2; if diag_index ==1 for i=1:n Z(i,:) = Z(i,:)/w(i); end for j=1:n Z(:,j) = Z(:,j)/w(j); end else Z = W_inv *Z* W_inv; end H = Z - H; i = 1; while (i<=k) %%% AA^* is not the identity matrix Ax(i) = H(I(i),J(i)); Ax(i) = Ax(i) + 1.0e-10*x(i); i = i+1; end end end end return %%% End of Jacobian_matrix.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% PCG method %%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% This is exactly the algorithm given by Hestenes and Stiefel (1952) %% An iterative method to solve A(x) =b %% The symmetric positive definite matrix M is a preconditioner for A. %% See Pages 527 and 534 of Golub and va Loan (1996) function [p,flag,relres,iterk] = pre_cg(b,I,J,tol,maxit,Omega12,WP,c,diag_index,W_inv,w_vec) n = length(WP); k1 = length(b); % initial x0=0 r = b; n2b = norm(b); % norm of b tolb = max(tol,min(0.1,n2b))*n2b; % relative tolerance if n2b > 1.0e2 maxit = min(1,maxit); end flag = 1; iterk = 0; relres = 1000; % to give a big value on relres p = zeros(k1,1); z = r./c; % z = M\r; here M =diag(c); if M is not the identity matrix rz1 = r'*z; rz2 = 1; d = z; for k = 1:maxit if k > 1 beta = rz1/rz2; d = z + beta*d; end w = Jacobian_matrix(d,I,J,Omega12,WP,diag_index,W_inv,w_vec); % w = A(d) if k1>n w = w + 1.0e-2*min(1.0, 0.1*n2b)*d; % perturb it to avoid numerical singularity end denom = d'*w; iterk = k; relres = norm(r)/n2b; % relative residue = norm(r)/norm(b) if denom <= 0 % sssss = 0 p = d/norm(d); % d is not a descent direction break else alpha = rz1/denom; p = p + alpha*d; r = r - alpha*w; end z = r./c; if norm(r) <= tolb iterk = k; relres = norm(r)/n2b; % relative residue = norm(r)/norm(b) flag =0; break end rz2 = rz1; rz1 = r'*z; end return %%% End of pre_cg.m %%% To generate the (approximate) diagonal preconditioner function c = precond_matrix(I,J,Omega12,WP) n = length(WP); k = length(I); [r,s] = size(Omega12); c = ones(k,1); H = WP'; H = H.*H; const_prec = 1; if (r0) if (k<=const_prec*n) % compute the exact diagonal preconditioner Ind = find(I~=J); k1 = length(Ind); if (k1>0) H1 = zeros(n,k1); for i=1:k1 H1(:,i) = WP(I(Ind(i)),:)'.*WP(J(Ind(i)),:)'; end end if (r0) H12_1 = H1(1:r,:)'*Omega12; end d = ones(r,1); j=0; for i=1:k if (I(i)==J(i)) c(i) = sum(H(1:r,I(i)))*(d'*H(1:r,J(i))); c(i) = c(i) + 2.0*(H12(I(i),:)*H(r+1:n,J(i))); else j=j+1; c(i) = sum(H(1:r,I(i)))*(d'*H(1:r,J(i))); c(i) = c(i) + 2.0*(H12(I(i),:)*H(r+1:n,J(i))); c(i) = c(i) + sum(H1(1:r,j))*(d'*H1(1:r,j)); c(i) = c(i) + 2.0*(H12_1(j,:)*H1(r+1:n,j)); c(i) = 0.5*c(i); end if c(i) < 1.0e-8 c(i) = 1.0e-8; end end else % if r>=n/2, use a complementary formula Omega12 = ones(r,s)-Omega12; H12 = Omega12*H(r+1:n,:); if(k1>0) H12_1 = Omega12*H1(r+1:n,:); end d = ones(s,1); dd = ones(n,1); j=0; for i=1:k if (I(i)==J(i)) c(i) = sum(H(r+1:n,I(i)))*(d'*H(r+1:n,J(i))); c(i) = c(i) + 2.0*(H(1:r,I(i))'*H12(:,J(i))); alpha = sum(H(:,I(i))); c(i) = alpha*(H(:,J(i))'*dd) - c(i); else j=j+1; c(i) = sum(H(r+1:n,I(i)))*(d'*H(r+1:n,J(i))); c(i) = c(i) + 2.0*(H(1:r,I(i))'*H12(:,J(i))); alpha = sum(H(:,I(i))); c(i) = alpha*(H(:,J(i))'*dd) - c(i); tmp = sum(H1(r+1:n,j))*(d'*H1(r+1:n,j)); tmp = tmp + 2.0*(H1(1:r,j)'*H12_1(:,j)); alpha = sum(H1(:,j)); tmp = alpha*(H1(:,j)'*dd) - tmp; c(i) = (tmp + c(i))/2; end if c(i) < 1.0e-8 c(i) = 1.0e-8; end end end else % approximate the diagonal preconditioner HH1 = H(1:r,:); HH2 = H(r+1:n,:); if (r