Efficiently repeating nullspace operation

I need to repeat alot of null-space operations (of 3 x 4 matrices) , and is looking for a way to improve speed. I am trying to work with matrices, and essentially what I want to do, if possible, is to increase speed between "tic" and "toc" in the following:

n=1e6;
Trow1=rand(n,4);
Trow2=rand(n,4);
Trow3=rand(n,4);
tic
Vh=zeros(n,4);
for i=1:n
    Vh(i,:)=null([ Trow1(i,:); Trow2(i,:); Trow3(i,:)] );
end  
toc

I am not looking for methods involving parallelization, but was hoping that I somehow could avoid the for-loop, or similar.

Answer 1

[Bruno Luong]

Here is one more method, so far fastest and derived from a much more direct calculation.

The idea is like this : For a fix matrix A (4x3), Vh is null(A') is a vector of norm_l2=1 that maximizes det([A, Vh]), I assume A is fullrank though.

Cauchy–Schwarz's inequality proof tells us that it maximizes when Vh is correlated to 3x3 sub-determinants of A. So I use to compute directly Vh.

The idea can be extended for any dimension array p x (p-1). For example for p==3, we gets the formula of the cross-product.

n=1e5;
Trow1=rand(n,4);
Trow2=rand(n,4);
Trow3=rand(n,4);

tic
T1 = reshape(Trow1.',[4 1 n]);
T2 = reshape(Trow2.',[4 1 n]);
T3 = reshape(Trow3.',[4 1 n]);
T = cat(2,T1,T2,T3);
[p,q,n] = size(T);
rdet(q);
c = 1:q;
Vh = zeros(p,n);
for i=1:p
    ri = setdiff(1:p,i);
    Vh(i,:) = (-1)^(i-1)*rdet(T,ri,c);
end
Vh = (Vh ./ sqrt(sum(Vh.^2,1))).';
toc % Elapsed time is 0.099139 seconds. (PS: run on different computer than before)

% Compute determinant recursively with bookeeping
function d = rdet(A,r,c)
persistent RC D

if nargin == 1
    % Reset data-base
    q = A;
    RC = cell(1,q);
    for k=1:q
        RC{k} = zeros(0,2*k);
    end
    D = cell(1,q);
    return
end

q = length(r);
if q == 1
    d = A(r,c,:);
    return
end

[tf,loc] = ismember([r,c],RC{q},'rows');
if tf
    d = D{q}(loc,1,:);
    return
end

% Update the row-column database
RC{q}(end+1,:) = [r,c];

% Recursive call of determinant
A1 = A(r,c(1),:);
c(1) = [];
di = zeros(size(A1));
for i=1:q
    ri = r;
    ri(i) = [];
    di(i,1,:) = rdet(A,ri,c);
end
A1(2:2:end,1,:) = -A1(2:2:end,1,:);
d = sum(A1.*di,1);

% Store the result
D{q}(end+1,1,:) = d;

end

Answer 2

[John D'Errico]

First, don't store your data in three separate arrays. Learn to use a 3-d array instead. That avoids having to aggregate those rows into one array. When you are talking about such small arrays, the time required for the null call is significant, but you need to recognize just how much time you are tossing in the bit bucket by what may seem minor inefficiencies.

n=1e5;
Trow1=rand(n,4);
Trow2=rand(n,4);
Trow3=rand(n,4);
tic
Vh=zeros(n,4);
for i=1:n
    Vh(i,:)=null([ Trow1(i,:); Trow2(i,:); Trow3(i,:)] );
end  
toc
Elapsed time is 4.347934 seconds.

Hey, its an old computer here, so I dropped n to 1e5.

T = rand(3,4,n);
tic
Vh=zeros(n,4);
for i=1:n
    Vh(i,:)=null(T(:,:,i));
end
toc
Elapsed time is 3.586421 seconds.

So, just using a rational way to store your data saved a bit. But there are a variety of ways to compute the null space of a matrix. Null is nice, because it is stable due to the use of an SVD. It also handles matrices of any size and rank.

tic
Vh=zeros(n,4);
for i=1:n
    [Q,~,~] = qr(T(:,:,i)');
    Vh(i,:) = Q(:,end)';
end
toc
Elapsed time is 1.078097 seconds.

Using a pivoted QR should be numerically stable, and much more efficient since we know we always want that last column of Q. By the way, you need to call qr as I did, even though you want only Q, because qr adjusts its behavior based on the number of output arguments.

Be careful though. Your code can easily fail as you wrote it. If those three rows are ever linearly dependent, then the null space will have dimension 2. In that case, the null space will be ambiguous.

[Bruno Luong]

+1

The next optimization step would be write down some sort of vectorized qr algorithm that can apply on a stack of matrices, using for example successive Housholder transformation as the calculations follow exactly the same path (in the case of non-pivoting) for all matrices.

[John D'Errico]

Yes, once you start thinking about it, there are surely ways to squeeze a bit more out of this.

Answer 3

[Bruno Luong]

Faster method still, solve a linear system by introducing an arbitrary vectors.

It needs such multiple system solver such as this FEX

n=1e5;
Trow1=rand(n,4);
Trow2=rand(n,4);
Trow3=rand(n,4);

tic
Vh=zeros(n,4);
for i=1:n
    Vh(i,:)=null([ Trow1(i,:); Trow2(i,:); Trow3(i,:)] );
end  
toc % Elapsed time is 3.044425 seconds.

tic
T1 = reshape(Trow1.',[4 1 n]);
T2 = reshape(Trow2.',[4 1 n]);
T3 = reshape(Trow3.',[4 1 n]);
T = cat(2,T1,T2,T3);
Vh = zeros(4,n);
for k=1:n
    [Q,~] = qr(T(:,:,k));
    Vh(:,k) = Q(:,4);
end
Vh = Vh.';
toc % Elapsed time is 0.882924 seconds.

tic
T1 = reshape(Trow1.',[1 4 n]);
T2 = reshape(Trow2.',[1 4 n]);
T3 = reshape(Trow3.',[1 4 n]);
T4 = randn(1,4,n);
T = cat(1,T1,T2,T3,T4);
Vh = MultiSolver(T,[0;0;0;1]); % FEX
n2 = sum(Vh.^2,1);
Vh = (Vh./sqrt(n2)).'; % normalize and transpose
toc % Elapsed time is 0.438108 seconds.

Answer 4

[Bruno Luong]

Here we go, a vectorized Housholder QR. You might need to replace the auto-expension with bsxfun if you run old MATLAB version.

tic
  T1 = reshape(Trow1.',[4 1 n]);
  T2 = reshape(Trow2.',[4 1 n]);
  T3 = reshape(Trow3.',[4 1 n]);
  T4 = randn([4 1 n]);
  R = cat(2,T1,T2,T3,T4);
  Q = repmat(eye(4),[1 1 n]);
  for k=1:size(R,2)
      i = k:size(R,1);
      j = k:size(R,2);
      u = R(i,k,:);
      a = sqrt(sum(u.^2,1));
      u1 = u(1,:,:);
      b = u1>0;
      a(b) = -a(b);
      u(1,:,:) = u1-a;
      nu = sqrt(sum(u.^2,1));
      nu(nu==0) = 1;
      v = u ./ nu;
      w = -2*v;
      v = reshape(v,[], 1, n);
      Q(i,:,:) = Q(i,:,:) + w .* sum(Q(i,:,:).*v,1);
      R(i,j,:) = R(i,j,:) + w .* sum(R(i,j,:).*v,1);
  end
  Vh = permute(Q(end,:,:),[3 2 1]);
toc % Elapsed time is 0.373327 seconds.

[Bruno Luong]

For complex data the loop must be modified to

      [p,q,n] = size(R);
      for k=1:q
          i = k:p;
          j = k:q;
          u = R(i,k,:);
          a = sqrt(sum(u.*conj(u),1));
          u1 = u(1,:,:);
          theta = angle(u1);
          a = a .* (-exp(1i*theta));
          u(1,:,:) = u1-a;
          nu = sqrt(sum(u.*conj(u),1));
          nu(nu==0) = 1;
          v = u ./ nu;
          w = -2*v;
          v = reshape(conj(v),[], 1, n);
          Q(i,:,:) = Q(i,:,:) + w .* sum(Q(i,:,:).*v,1);
          R(i,j,:) = R(i,j,:) + w .* sum(R(i,j,:).*v,1);
      end
      Vh = permute(conj(Q(end,:,:)),[3 2 1])