There are many methods to obtain null vectors to a given set of vectors. One of the numerically best methods is implemented in MATLAB, based on singular value decomposition, SVD. Null vectors can be used for many different things; one of them is to remove properties that are ore orthogonal to another, given property. Such operations are used often within chemometrics and one of the algorithms that got most attention recently is the OPLS and the O2PLS algorithms by Trygg and Wold, but I wanted to use it together with my research on PLS-algorithms.
Example problems:
- You have a vector w and want to find all vectors that are orthogonal to it.
- You have a matrix W with orthonormal vectors and want to find all vectors that are orthogonal to them,
I wanted to have a method that is similar to the NULL-function that can be found in MATLAB, meaning that I wanted to have all the null vectors available. If I would need only one or two more null vectors, other methods could have been used. The ones that seemed to be most suitable were the ones based on Householder reflections.
The data that I wanted to pad null vectors to consist of numerically highly orthonormal vectors. I therefore thought I could make something faster than MATLAB’s method based on SVD. I tried many ways but one of the most promising methods was a modified QR-decomposition, using its full mode version. After getting looping right, it turned out that the QR decomposition implemented in MATLAB was unbeatable, even after porting it to C++ and converting it to MATLAB executables. I therefore gave up on making my own modified QR-decomposition and used the one implemented in MATLAB.
Next surprise came when I compared the execution times and the orthogonaly properties when using the null space vectors obtained from MATLAB’s SVD with the ones from MATLAB’s QR. The calculation time was very similar for the two. The most logical choise would perhaps then be to use the one that already existed, the one based on SVD in MATLAB. However, because the SVD-algorithm is complicated, I decided to use the QR version, which would also be much easier to implement in other code, for example C++.
After this work, I thought I would not reach any further in the null space and decided to contine the work on a new PLS algorithm instead.
My null space algorithm based on QR
%nullA = nullQR(A)
%
% get the null space of A by continuing using QR decomposition,
% *** assuming that: ***
% *** m > n ***
% *** A is already perfectly orthogonal. ***
% No check will be done so use with care.
% The I/O is different to NULL using a transposed matrix.
%
% See also NULL, SVD, ORTH, RANK, RREF.
%
%
% Martin Andersson, 2011-11
[m,n] = size(A);
[Q,R] = qr(A);
Z = Q(:,n+1:m);
Null space from MATLAB’s SVD
%NULL Null space.
% Z = NULL(A) is an orthonormal basis for the null space of A obtained
% from the singular value decomposition. That is, A*Z has negligible
% elements, size(Z,2) is the nullity of A, and Z’*Z = I.
%
% Z = NULL(A,‘r’) is a “rational” basis for the null space obtained
% from the reduced row echelon form. A*Z is zero, size(Z,2) is an
% estimate for the nullity of A, and, if A is a small matrix with
% integer elements, the elements of R are ratios of small integers.
%
% The orthonormal basis is preferable numerically, while the rational
% basis may be preferable pedagogically.
%
% Example:
%
% A =
%
% 1 2 3
% 1 2 3
% 1 2 3
%
% Z = null(A);
%
% Computing 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 support for input A:
% float: double, single
%
% See also SVD, ORTH, RANK, RREF.
% Copyright 1984–2006 The MathWorks, Inc.
m,n] = size(A);
if (nargin > 1) && isequal(how,‘r’)
% Rational basis
[R,pivcol] = rref(A);
r = length(pivcol);
nopiv = 1:n;
nopiv(pivcol) = [];
Z = zeros(n,n‑r,class(A));
if n > r
Z(nopiv,:) = eye(n‑r,n‑r,class(A));
if r > 0
Z(pivcol,:) = ‑R(1:r,nopiv);
end
end
else
% Orthonormal 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