Using safe_eig() code to prevent rare eigenvalue decomp. bug

private/eig_backup.m

@@ -0,0 +1,22 @@
+function [V,D] = eig_backup( X )
+% [V,D] = eig_backup(X)
+%   computes the eigenvalue decomposition of a symmetric matri
+%   using the SVD. This is designed for cases when the main eig()
+%   routine doesn't work due to the bug described at 
+%       http://ask.cvxr.com/t/eig-did-not-converge-in-prox-trace/996/4
+%
+% This uses the SVD, so a bit slower sometimes...
+% July 13 2015
+
+[V,D,W] = svd(X);
+
+
+d = diag(D).' .* sign(real(dot(V,W,1)));
+% and sort it to conform to Matlab's order...
+[d,ind]     = sort(d);
+D           = diag(d);
+V           = V(:,ind);
+
+% Another option, but typically slow...
+%   (calls generalized eigenvalue decomp.)/
+% [V,D]   = eig(X,eye(size(X)));

private/safe_eig.m

@@ -0,0 +1,14 @@
+function [V,D] = safe_eig(X)
+% [V,D] = safe_eig(X)
+%   calls [V,D] = eig(X) and tests for a common error
+%   (http://ask.cvxr.com/t/eig-did-not-converge-in-prox-trace/996/4)
+%   and if it finds it, runs a replacement algorithm
+try
+    [V,D]   = eig(X);
+catch ME
+    if (strcmpi(ME.identifier,'MATLAB:eig:NoConvergence'))
+        [V,D]   = eig_backup(X);
+    else
+        rethrow(ME);
+    end
+end

proj_maxEig.m

@@ -45,7 +45,7 @@
 v = 0;
 if nargin > 3 && t > 0,
     if SP, X = full(X); end % svd, eig, and norm require full matrices
-    [V,D]   = eig(X); % just in case X is sparse
+    [V,D]   = safe_eig(X); % just in case X is sparse
     [dum,D] = eproj(diag(D),q); % not q*t, since we are just projecting...
     X       = V*diag(D)*V';
     if SP, X = sparse(X); end
@@ -84,7 +84,7 @@
     while ~ok
         K = min( [K,M,N] );
         if K > min(M,N)/2 || K > (min(M,N)-2) || min(M,N) < 20
-            [V,D]   = eig(full((X+X')/2));
+            [V,D]   = safe_eig(full((X+X')/2));
             ok = true;
         else
             [V,D] = eigs( X, K, SIGMA, opts );

proj_psd.m

@@ -37,7 +37,7 @@
 	v = 0;
     X = full(X+X'); % divide by 2 later
     if isReal, X = real(X); end
-    [V,D]=eig(X); % we don't yet take advantage of sparsity here
+    [V,D]=safe_eig(X); % we don't yet take advantage of sparsity here
     D  = max(0.5*diag(D),0);
     tt = D > 0;
     V  = bsxfun(@times,V(:,tt),sqrt(D(tt,:))');
@@ -84,7 +84,7 @@
     while ~ok
         K = min( [K,N] );
         if K > N/2 || K > N-2 || N < 20
-            [V,D]   = eig(full((X+X')/2));
+            [V,D]   = safe_eig(full((X+X')/2));
             ok = true;
         else
             [V,D] = eigs( X, K, SIGMA, opts );

proj_psdUTrace.m

@@ -114,7 +114,7 @@
         ctr   = ctr + 1;
         K     = min( K, N );
         if K > N/2 || K > maxK
-            [V,D] = eig(X); ok = true;
+            [V,D] = safe_eig(X); ok = true;
             D     = diag(D);
             break;
         end
@@ -195,7 +195,7 @@
 end
 if nargin > 4 && t > 0,
     nCalls = nCalls + 1;
-    [V,D]=eig(X);
+    [V,D]=safe_eig(X);
     [dum,D] = eproj(diag(D),t);
     tt = D > 0;
     V  = bsxfun(@times,V(:,tt),sqrt(D(tt,:))');

proj_spectral.m

@@ -84,7 +84,7 @@
 v = 0;
 if nargin > 3 && t > 0,
     if SP, X = full(X); end % svd, eig, and norm require full matrices
-    [V,D]   = eig(X); % just in case X is sparse
+    [V,D]   = safe_eig(X); % just in case X is sparse
     [dum,D] = eproj(diag(D),q); % not q*t, since we are just projecting...
     X       = V*diag(D)*V';
     if SP, X = sparse(X); end
@@ -135,7 +135,7 @@
             opts.tol = 1e-3;
         end
         if K > min(M,N)/2
-            [V,D]   = eig(full((X+X')/2));
+            [V,D]   = safe_eig(full((X+X')/2));
             ok = true;
         end
     end

prox_maxEig.m

@@ -42,7 +42,7 @@
 X = (X+X')/2; % Matlab will make it exactly symmetric
 
 if nargin == 3 && t > 0,
-    [V,S]       = eig(full(X));
+    [V,S]       = safe_eig(full(X));
 %     op          = prox_linf(q); % for prox_spectral
     op          = prox_max(q);
     s           = diag(S);

prox_spectral.m

@@ -81,7 +81,7 @@
 X = (X+X')/2; % Matlab will make it exactly symmetric
 
 if nargin == 3 && t > 0,
-    [V,S]       = eig(full(X));
+    [V,S]       = safe_eig(full(X));
     op          = prox_linf(q);
     s           = diag(S);
     [dummy,s]   = op(s,t);

prox_trace.m

@@ -87,25 +87,9 @@
         if isReal
             % July 10 2015, adding try-catch statement
             % to deal with bug in eig() sometimes
-            try
-                [V,D]   = eig(real(full((X+X')/2)));
-            catch ME
-                if (strcmpi(ME.identifier,'MATLAB:eig:NoConvergence'))
-                    [V,D]   = eig(real(full((X+X')/2)),eye(size(X)));
-                else
-                    rethrow(ME);
-                end
-            end
+            [V,D]   = safe_eig(real(full((X+X')/2)));
         else
-            try
-                [V,D]   = eig(full((X+X')/2));
-            catch ME
-                if (strcmpi(ME.identifier,'MATLAB:eig:NoConvergence'))
-                    [V,D]   = eig(full((X+X')/2),eye(size(X)));
-                else
-                    rethrow(ME);
-                end
-            end
+            [V,D]   = safe_eig(full((X+X')/2));
         end
     else
 
@@ -131,15 +115,7 @@
             K = min( [K,M,N] );
 
             if K > min(M,N)/2
-                try
-                    [V,D]   = eig(full((X+X')/2));
-                catch ME
-                    if (strcmpi(ME.identifier,'MATLAB:eig:NoConvergence'))
-                        [V,D]   = eig(full((X+X')/2),eye(size(X)));
-                    else
-                        rethrow(ME);
-                    end
-                end
+                [V,D]   = safe_eig(full((X+X')/2));
                 ok = true;
                 break;
             end

smooth_logdet.m

@@ -63,7 +63,7 @@
         % the function is being used in a "nonsmooth" fashion
         % i.e. return g = argmin_g  -q*log(det(g)) + 1/(2t)||g-x||^2
         x = full(x+x')/2;  % March 2015, project it to be symmetric
-        [V,D]   = eig(x);
+        [V,D]   = safe_eig(x);
         d       = diag(D);
         % This is OK: input need not be pos def
         %if any(d<=0),
@@ -108,7 +108,7 @@
         % i.e. return g = argmin_g  -q*log(det(g)) + 1/(2t)||g-x||^2
         x       = x - t*C;  
         x = full(x+x')/2;  % March 2015, project it to be symmetric
-        [V,D]   = eig(x);
+        [V,D]   = safe_eig(x);
         d       = diag(D);
         % This is OK: input need not be pos def
         l       = ( d + sqrt( d.^2 + 4*t*q ) )/2;

smooth_quad.m

@@ -73,7 +73,7 @@
     else
         P = 0.5 * ( P + P' );
         if use_eig
-            [V,DD] = eig(P); dd = diag(DD);
+            [V,DD] = safe_eig(P); dd = diag(DD);
         else
             V = []; dd = [];
         end