Numerically good null spaces

There are many meth­ods to obtain null vec­tors to a given set of vec­tors. One of the numer­i­cally best meth­ods is imple­mented in MATLAB, based on sin­gu­lar value decom­po­si­tion, SVD. Null vec­tors can be used for many dif­fer­ent things; one of them is to remove prop­er­ties that are ore orthog­o­nal to another, given prop­erty. Such oper­a­tions are used often within chemo­met­rics and one of the algo­rithms that got most atten­tion recently is the OPLS and the O2PLS algo­rithms by Trygg and Wold, but I wanted to use it together with my research on PLS-algorithms.

Exam­ple problems:

  • You have a vec­tor w and want to find all vec­tors that are orthog­o­nal to it.
  • You have a matrix W with ortho­nor­mal vec­tors and want to find all vec­tors that are orthog­o­nal to them,

I wanted to have a method that is sim­i­lar to the NULL-func­tion that can be found in MATLAB, mean­ing that I wanted to have all the null vec­tors avail­able. If I would need only one or two more null vec­tors, other meth­ods could have been used. The ones that seemed to be most suit­able were the ones based on House­holder reflections.

The data that I wanted to pad null vec­tors to con­sist of numer­i­cally highly ortho­nor­mal vec­tors. I there­fore thought I could make some­thing faster than MATLAB’s method based on SVD. I tried many ways but one of the most promis­ing meth­ods was a mod­i­fied QR-decom­po­si­tion, using its full mode ver­sion. After get­ting loop­ing right, it turned out that the QR decom­po­si­tion imple­mented in MATLAB was unbeat­able, even after port­ing it to C++ and con­vert­ing it to MATLAB exe­cuta­bles. I there­fore gave up on mak­ing my own mod­i­fied QR-decom­po­si­tion and used the one imple­mented in MATLAB.

Next sur­prise came when I com­pared the exe­cu­tion times and the orthog­o­naly prop­er­ties when using the null space vec­tors obtained from MATLAB’s SVD with the ones from MATLAB’s QR. The cal­cu­la­tion time was very sim­i­lar for the two. The most log­i­cal choise would per­haps then be to use the one that already existed, the one based on SVD in MATLAB. How­ever, because the SVD-algo­rithm is com­pli­cated, I decided to use the QR ver­sion, which would also be much eas­ier to imple­ment in other code, for exam­ple C++.

After this work, I thought I would not reach any fur­ther in the null space and decided to con­tine the work on a new PLS algo­rithm instead.

My null space algo­rithm based on QR

func­tion Z = nul­lQR(A);
%nullA = nullQR(A)
%
% get the null space of A by con­tin­u­ing using QR decomposition,
% *** assum­ing that: ***
% ***     m > n ***
% ***     A is already per­fectly orthogonal. ***
% No check will be done so use with care.
% The I/O is dif­fer­ent to NULL using a trans­posed matrix.
%
%   See also NULL, SVD, ORTH, RANK, RREF.
%
%
% Mar­tin Ander­s­son, 2011-11
[m,n] = size(A);
[Q,R] = qr(A);
Z = Q(:,n+1:m);

Null space from MAT­LAB’s SVD

func­tion Z = null(A,how)
%NULL   Null space.
%   Z = NULL(A) is an ortho­nor­mal basis for the null space of A obtained
%   from the sin­gu­lar value decom­po­si­tion.  That is,  A*Z has negligible
%   ele­ments, size(Z,2) is the nul­lity of A, and Z’*Z = I.
%
%   Z = NULL(A,‘r’) is a “ratio­nal” basis for the null space obtained
%   from the reduced row ech­e­lon form.  A*Z is zero, size(Z,2) is an
%   esti­mate for the nul­lity of A, and, if A is a small matrix with
%   inte­ger ele­ments, the ele­ments of R are ratios of small integers.
%
%   The ortho­nor­mal basis is prefer­able numer­i­cally, while the rational
%   basis may be prefer­able pedagogically.
%
%   Exam­ple:
%
%       A =
%
%           1     2     3
%           1     2     3
%           1     2     3
%
%       Z = null(A);
%
%       Com­put­ing the 1‑norm of the matrix A*Z will be
%       within a small tolerance
%
%       norm(A*Z,1)< 1e-12
%       ans =
%
%          1
%
%       null(A,‘r’) =
%
%          ‑2    ‑3
%           1     0
%           0     1
%
%   Class sup­port for input A:
%      float: dou­ble, single
%
%   See also SVD, ORTH, RANK, RREF.
%   Copy­right 1984–2006 The Math­Works, Inc.
m,n] = size(A);
if (nar­gin > 1) && ise­qual(how,‘r’)
% Ratio­nal basis
[R,pivcol] = rref(A);
r = length(piv­col);
nopiv = 1:n;
nopiv(piv­col) = [];
Z = zeros(n,n‑r,class(A));
if n > r
Z(nopiv,:) = eye(n‑r,n‑r,class(A));
if r > 0
Z(piv­col,:) = ‑R(1:r,nopiv);
end
end
else
% Ortho­nor­mal basis
[~,S,V] = svd(A,0);
if m > 1, s = diag(S);
elseif m == 1, s = S(1);
else s = 0;
end
tol = max(m,n) * max(s) * eps(class(A));
r = sum(s > tol);
Z = V(:,r+1:n);
end

Leave a Comment

Your email address will not be published. Required fields are marked *

This website uses cookies. By continuing to use this site, you accept our use of cookies. 

Scroll to Top