Vectorizing code for sign flipping of eigenvectors

So I have this function to resolve the sign ambiguity of eigenvector decomposition (i.e. to remove side effects from eig). Since the code still includes for loops, it's running really slow. I'm wondering if someone could give my any tips on how to vectorize this following function. I'd be grateful, even if it's only part of it.

function [v_flip, w_flip] = sign_flip(x,v,d,w)
    v_flip = zeros(size(v));
    w_flip = zeros(size(w));

    [I, J] = size(x);
    
    for k = 1:size(d,1)
       tmp_sum = zeros(I,J);
       
       for m = 1:length(d)
          if m ~= k
           tmp_sum = tmp_sum + d(m,m)*v(:,m)*w(:,m)';
          end
       end

       y = x - tmp_sum;
        
       sk_left = sum(sign(v(:,k)'*y).*(v(:,k)'*y).^2);
       sk_right = sum(sign(w(:,k)'*y').*(w(:,k)'*y').^2);
        
       if sk_left*sk_right < 0
          if sk_left < sk_right
              sk_left = -sk_left;
          else
              sk_right = -sk_right;
          end
       end
       
       v_flip(:,k) = sign(sk_left)*v(:,k);
       w_flip(:,k) = sign(sk_right)*w(:,k);
    end
end

And a demonstation of how the function works:

A = randi(5,[5,5]);
A = A'*A;
[v,d,w] = eig(A);
v =
   -0.1038    0.2939   -0.6668    0.4562    0.5001
   -0.4406   -0.0441    0.6623    0.4226    0.4321
   -0.2923    0.3193    0.0117   -0.7696    0.4693
   -0.0042   -0.8919   -0.2125   -0.1448    0.3720
    0.8424    0.1194    0.2672    0.0095    0.4523

[v_flip, w_flip] = sign_flip(A,v,d,w);
v_flip =
   -0.1038   -0.2939   -0.6668   -0.4562    0.5001
   -0.4406    0.0441    0.6623   -0.4226    0.4321
   -0.2923   -0.3193    0.0117    0.7696    0.4693
   -0.0042    0.8919   -0.2125    0.1448    0.3720
    0.8424   -0.1194    0.2672   -0.0095    0.4523

Answer 1

[John D'Errico]

I would suggest your code is needlessly complex. Trivial to write too.

Choose an orientation for the eigenvectors such that the largest element in absolute value always has a positive sign. (Or something like that. You might also choose to fix the sign to be positive of the first element in the vector that is distinct from zero by some tolerance. )

A = rand(5); A = A'*A;

[V,D] = eig(A);
tol = 1e-15;
t = abs(V) > tol;
[~,ind] = max(t,[],1);
% ind will usually be just a vector of 1's
% the rare case when it is not is when the
% first element of an eigenvector is smaller
% than tol in absolute value.
n = size(A,1);
S = sign(V(ind + (0:n:n^2-1)));

% Now just flip the signs. I'll use a trick that
% allowed in R2016b or later. So if you have an
% older release, you would use bsxfun.
Vflip = V.*S;

Did it work? Of course.

S
S =
     1     1    -1     1     1     
     
V
V =
      0.59868      0.38023     -0.39213       0.1419      0.56842
      0.08037     -0.47118      0.49184      0.59306      0.42179
     -0.76384      0.33681     -0.26517      0.37595      0.30241
     -0.11257      0.42384       0.6502     -0.47288      0.40164
     -0.19749     -0.58338     -0.33354     -0.51302      0.49621
     
Vflip
Vflip =
      0.59868      0.38023      0.39213       0.1419      0.56842
      0.08037     -0.47118     -0.49184      0.59306      0.42179
     -0.76384      0.33681      0.26517      0.37595      0.30241
     -0.11257      0.42384      -0.6502     -0.47288      0.40164
     -0.19749     -0.58338      0.33354     -0.51302      0.49621