EE5120 Linear Algebra - Chapter 3 Matlab resources

Companion page: Chapter 3 — Orthogonality — the notes that this page goes with.

Commands for the material of Chapter 3: norms, projections, least squares, orthogonal matrices and \(QR\). Everything works in GNU Octave too, and the places where the two differ are called out. The running example is the velocity fit of the notes,

t = [-1; 1; 2];
A = [ones(3,1) t];      % 3 x 2, independent columns
b = [1; 1; 3];

Norms

x = [1; 0; 2; 2];
norm(x)                  % 3, the Euclidean norm: sqrt(1+0+4+4)
norm(x,2)                % same thing
norm(x,1)                % 5,  sum of absolute values
norm(x,Inf)              % 2,  largest absolute value
nnz(x)                   % 3,  the "L0 count" -- note it is not called a norm

norm on a matrix is not the Frobenius norm by default — it is the largest amount by which \(A\) can stretch a vector, which needs the SVD:

norm(A)                  % the 2-norm (largest singular value)
norm(A,'fro')            % the Frobenius norm of the notes
sqrt(sum(A(:).^2))       % same thing, computed by hand

Watch the homogeneity failure of the \(\ell_0\) count, which is what disqualifies it as a norm:

[nnz(x), nnz(2*x)]       % 3 3 -- unchanged, whereas a norm would double
[norm(x), norm(2*x)]     % 3 6 -- a real norm doubles

Cauchy—Schwarz and the angle between two vectors:

u = [1; 2; 3]; v = [4; 5; 6];
abs(u'*v) <= norm(u)*norm(v)          % 1, always
acos( u'*v/(norm(u)*norm(v)) )        % the angle in radians, 0.2257
subspace(u,v)                         % same angle, computed stably

Orthogonal complements and the splitting

null and orth return orthonormal bases, which is exactly what this chapter is about.

M = [1 3 3 2; 2 6 9 7; -1 -3 3 4];    % the Chapter 2 example, rank 2
Qr = orth(M');                        % orthonormal basis of the row space
Nu = null(M);                         % orthonormal basis of the null space
[size(Qr,2) size(Nu,2)]               % 2 and 2, adding to n = 4
norm(Qr'*Nu)                          % ~1e-16: the two spaces are orthogonal
rank([Qr Nu])                         % 4 -- together they span all of R^4

Now split a vector in two, and check the pieces are unique and Pythagorean:

x  = [3; -1; 4; 2];
xr = Qr*(Qr'*x);          % row space part
xn = Nu*(Nu'*x);          % null space part
norm(x - (xr+xn))         % ~1e-15: the two pieces add back to x
xr'*xn                    % ~1e-15: they are perpendicular
[norm(x)^2, norm(xr)^2 + norm(xn)^2]   % equal -- Pythagoras
norm(M*x - M*xr)          % ~1e-14: M discards the null space part

Projection and least squares

Three ways to the same answer. Use the last one in real work:

xhat1 = (A'*A)\(A'*b)     % the normal equations, solved by elimination
xhat2 = pinv(A)*b         % via the pseudo-inverse
xhat3 = A\b               % backslash -- for a tall A this IS least squares
[xhat1 xhat2 xhat3]       % all three columns are (9/7, 4/7) = (1.2857, 0.5714)

Backslash on a tall system does not try to solve \(Ax=b\) — it minimises \(\|b-Ax\|\), and it does it by \(QR\) rather than by forming \(A^TA\). Never write inv(A'*A)*A'*b; it is slower and less accurate than all three of the above.

The projection, the error, and the two checks from the notes:

p = A*xhat1               % (5/7, 13/7, 17/7) = (0.7143, 1.8571, 2.4286)
e = b - p                 % (2/7, -6/7, 4/7)
A'*e                      % ~1e-15: e is perpendicular to both columns
norm(e)^2                 % 8/7 = 1.1429

The projection matrix, and its two defining properties:

P = A*inv(A'*A)*A';       % explicit inverse only because we want the matrix
norm(P*b - p)             % ~1e-16
norm(P*P - P)             % ~1e-16: idempotent
norm(P' - P)              % ~1e-16: symmetric
rank(P)                   % 2 = n, so P is NOT invertible on R^3
norm((eye(3)-P)*b - e)    % ~1e-16: I-P projects onto the complement

Projection onto a single line is the \(n=1\) case, and a*a' is a rank one matrix:

a  = [1; 1; 1];
Pa = (a*a')/(a'*a);       % rank 1
Pa*b                      % (5/3, 5/3, 5/3): the average, repeated
rank(Pa)                  % 1
norm(Pa*b) <= norm(b)     % 1 -- a projection never lengthens

Plot the fit and see both pictures of the notes:

tt = linspace(-1.5, 2.5, 100);
plot(t, b, 'ko', 'MarkerFaceColor','k'); hold on; grid on
plot(tt, xhat1(1) + xhat1(2)*tt, 'b-', 'LineWidth', 1.5);
for k = 1:3
  plot([t(k) t(k)], [b(k) p(k)], 'r--', 'LineWidth', 1.5);   % the errors
end
xlabel('t'); ylabel('v'); legend('data','least squares fit','error')

Weighted least squares

Distrust the third reading by a factor of ten in variance:

w = [1; 1; 0.1];                  % weights, 1/variance
W = diag(sqrt(w));
xw = (W*A)\(W*b)                  % (1.0800, 0.1600) = (27/25, 4/25)
xwn = (A'*diag(w)*A)\(A'*diag(w)*b);   % the normal equations version
norm(xw - xwn)                    % ~1e-16: same answer

The unweighted slope was \(4/7 = 0.5714\); distrusting the third reading, which is the one high point, drops it to \(4/25 = 0.16\). Both lscov(A,b,w) and the two lines above give the same thing, in Matlab and in Octave alike.

Under-determined systems

Here backslash and pinv genuinely disagree, and the difference is the point of the section.

F = [1 2 3; 4 5 7];       % 2 x 3, fat, independent rows
c = [6; 15];
x_bs  = F\c
x_min = pinv(F)*c         % (3/7, 6/7, 9/7) = (0.4286, 0.8571, 1.2857)
x_frm = F'*((F*F')\c)     % the formula of the notes -- same as pinv
norm(x_min - x_frm)       % ~1e-16
norm(F*x_bs - c)          % ~1e-15: both are exact solutions ...
norm(F*x_min - c)         % ~1e-15: ... the error is zero either way
[norm(x_bs) norm(x_min)]

Important

This is where the two part company, and neither is wrong.

Matlab returns a basic solution — \((0.6,\; 0,\; 1.8)\), with as many zeros as it can manage — of length \(1.8974\). It prints no warning, because the matrix is not rank deficient.

Octave returns the minimum norm solution, silently: the same \((3/7,\; 6/7,\; 9/7)\) that pinv gives, of length \(1.6036\).

pinv agrees on both. When you want a particular solution, ask for it by name rather than trusting backslash.

Confirm that the minimum norm solution is the one lying in the row space, i.e. the one with no null space component:

norm(null(F)'*x_min)      % ~1e-16 on both: nothing in the null space
norm(null(F)'*x_bs)       % 1.0142 in Matlab; ~1e-16 in Octave

Orthogonal matrices

th = pi/6;
Q  = [cos(th) -sin(th); sin(th) cos(th)];
norm(Q'*Q - eye(2))       % ~1e-16
norm(Q*Q' - eye(2))       % ~1e-16 too, because Q is square
x  = [3; 4];
[norm(x) norm(Q*x)]       % both 5: length is preserved

Permutation matrices are orthogonal, which is the Chapter 1 fact reappearing:

Pm = eye(4); Pm = Pm([3 1 4 2],:);
norm(Pm'*Pm - eye(4))     % 0, exactly -- no round off at all here

A tall \(Q\) has a left inverse and no right inverse:

Qt = orth(A);             % 3 x 2, orthonormal columns
norm(Qt'*Qt - eye(2))     % ~1e-16
norm(Qt*Qt' - eye(3))     % about 1 -- NOT the identity
rank(Qt*Qt')              % 2, so it cannot be I(3)
norm(Qt*Qt' - P)          % ~1e-15: it is the same projection matrix as above
norm(Qt*(Qt'*b) - p)      % ~1e-15: least squares by inner products alone

Gram—Schmidt and \(QR\)

G = [1 1 2; 0 0 1; 1 0 0];        % the worked example of the notes
[Q,R] = qr(G);
norm(Q*R - G)                     % ~1e-15
norm(Q'*Q - eye(3))               % ~1e-16
norm(triu(R) - R)                 % 0: R is upper triangular

Important

The \(Q\) you get back will not match the one in the notes entry by entry. In the versions tested (Matlab R2025b and Octave 9.1) the first two columns come out negated; the sign convention is an implementation detail and can differ between releases, so test the factorisation rather than the entries. \(QR\) is only unique once you fix the signs of the diagonal of \(R\), and library routines use Householder reflections rather than Gram—Schmidt, so they make a different choice. This is the same situation as null in Chapter 2: test the factorisation, not the entries.

d  = sign(diag(R));               % the notes take all of these positive
Qn = Q*diag(d); Rn = diag(d)*R;   % renormalised: now it matches the notes
norm(Qn*Rn - G)                   % ~1e-15, still a valid factorisation
Rn                                % diagonal is sqrt(2), 1/sqrt(2), 1 as in the notes

Gram—Schmidt written out, to see that it is the process of the notes and not magic:

function [Q,R] = mygs(A)          % classical Gram-Schmidt
  [m,n] = size(A); Q = zeros(m,n); R = zeros(n,n);
  for j = 1:n
    v = A(:,j);
    for i = 1:j-1
      R(i,j) = Q(:,i)'*A(:,j);
      v = v - R(i,j)*Q(:,i);
    end
    R(j,j) = norm(v);
    Q(:,j) = v/R(j,j);
  end
end

[Qg,Rg] = mygs(G);
Qg                                % now this DOES match the notes
norm(Qg*Rg - G)                   % ~1e-16

Save the function in mygs.m, or paste it at the end of a script file in Matlab; Octave allows it anywhere.

Least squares by \(QR\), which is what backslash does internally:

[Qa,Ra] = qr(A,0);                % economy size: Qa is 3x2, Ra is 2x2
xq = Ra\(Qa'*b)                   % (9/7, 4/7) again
norm(xq - xhat1)                  % ~1e-15

Why bother, when the normal equations are shorter to write? Because squaring the matrix squares the conditioning:

cond(A), cond(A'*A)               % 1.8708 and 3.5000 -- exactly the square
H  = hilb(12);                    % the Hilbert matrix of the notes
cond(H)                           % ~1e16: numerically singular in double
[cond(hilb(3)) cond(hilb(6)) cond(hilb(10))]   % 524, 1.5e7, 1.6e13

That last line is the reason the monomials \(1, z, z^2, \ldots\) are a bad basis.

Function spaces

Inner products become integrals, and everything else is unchanged.

f = @(x) sin(x); g = @(x) cos(x);
integral(@(x) f(x).*g(x), 0, 2*pi)      % ~1e-16: orthogonal
integral(@(x) f(x).^2,    0, 2*pi)      % pi = 3.1416, so NOT orthonormal
integral(@(x) f(x).*sin(2*x), 0, 2*pi)  % ~1e-16: different harmonics too

Octave has both integral and quad, so either call works there. quad is the older name; Matlab still ships it, but it is documented as not recommended, so prefer integral in new code.

The Gram matrix of the monomials on \((0,1)\) is the Hilbert matrix, and its pivots are the Gram—Schmidt lengths of the notes:

H3 = hilb(3)                      % [1 1/2 1/3; 1/2 1/3 1/4; 1/3 1/4 1/5]
Rc = chol(H3);                    % H3 = Rc'*Rc, so Rc is the R of the notes
diag(Rc)'.^2                      % 1, 0.083333, 0.0055556 = 1, 1/12, 1/180
1./diag(Rc)'.^2                   % 1, 12, 180 -- easier to recognise

Do not reach for lu here. lu pivots rows, and on hilb(3) it exchanges two of them, so diag(U) comes back as \(1,\; 1/12,\; -1/180\) — the magnitudes survive but the sign is an artefact of the exchange. chol needs no pivoting on a matrix of this kind and returns the pivots directly as \(R_{jj}^2\), which is the theorem at the end of the notes.

Check that against Gram—Schmidt applied to \(\{1, z, z^2\}\) by sampling the interval finely, so that the integral becomes a dot product:

z  = linspace(0,1,20001)';  h = z(2)-z(1);
B  = [ones(size(z)) z z.^2];      % the three basis functions, sampled
Bs = B*sqrt(h);                   % scaling so that Bs'*Bs approximates the integrals
Bs'*Bs                            % the Hilbert matrix, to about four digits
[Qf,Rf] = mygs(Bs);
diag(Rf)'.^2                      % 1.00005, 0.083346, 0.0055569 -- the pivots again

The columns of Qf are the shifted Legendre polynomials, up to the sampling scale factor:

q1 = Qf(:,2)/sqrt(h);  q2 = Qf(:,3)/sqrt(h);
[q1(1) 2*sqrt(3)*(z(1)-0.5)]      % -1.7319 and -1.7321
max(abs(q2 - sqrt(5)*(6*z.^2 - 6*z + 1)))   % ~4e-4, the sampling error