Companion page: Chapter 2 — Vector spaces — the notes that this page goes with.

Commands for the material of Chapter 2: echelon form, the four subspaces, rank, and transformations. As before, everything works in GNU Octave too. The matrix used throughout is the running example of the notes,

A = [1 3 3 2; 2 6 9 7; -1 -3 3 4];   % 3 x 4, rank 2
b = [1; 5; 5];

Echelon form and rank

R = rref(A)              % reduced row echelon form
[R,p] = rref(A)          % p lists the pivot columns -- here [1 3]
rank(A)                  % 2
size(A,2) - rank(A)      % number of free variables, n - r = 2

There is no built-in for the hand echelon form \(U\) of the notes, and lu is not it — lu pivots rows, and on a rank-deficient matrix the \(U\) it returns need not be in echelon form at all. This matrix is exactly such a case:

[L,U,P] = lu(A)          % a correct factorisation, but NOT the U of the notes
P*A - L*U                % zero to round off: the factorisation itself is fine

You get \(U = \begin{bmatrix} 2 & 6 & 9 & 7 \\ 0 & 0 & -1.5 & -1.5 \\ 0 & 0 & 7.5 & 7.5\end{bmatrix}\), whose second and third rows have their leading entries in the same column, and which has no zero row. So use rref when you want the echelon story, and lu only when you want a factorisation.

Note rank does not count pivots by elimination — it counts singular values above a tolerance, which is the numerically sound way and is why it can disagree with a hand count on a badly scaled matrix. On this matrix the tolerance is irrelevant (rank(A,1e-12) is still 2); to see it bite, try rank([1 1; 1 1+1e-13]), which is 2 by default and 1 with a tolerance of 1e-12.

The null space

null(A)                  % orthonormal basis, 4 x 2 -- not the special solutions
null(A,'r')              % rational basis: these ARE the special solutions (Matlab only)
A*null(A)                % zero to round off

The two special solutions in the notes are \((-3,1,0,0)\) and \((1,0,-1,1)\). null(A,'r') returns exactly these; plain null(A) returns a different basis for the same subspace, orthonormal instead of rational. Both are correct — a basis is not unique. Octave has no 'r' option, so build them by hand from R there.

The column space

The pivot columns of \(A\) — not of R — form a basis:

[R,p] = rref(A);
basisC = A(:,p)          % columns 1 and 3 of A
orth(A)                  % an orthonormal basis for the same space
rank([A(:,p) A])         % still 2: the other columns add nothing

Confirm that elimination changes the column space while leaving the null space alone:

rank([A(:,1) R(:,1)])    % 2 -- so col 1 of A and of R span different lines
norm(A*null(rref(A)))    % ~1e-15: R's null vectors are also A's
norm(rref(A)*null(A))    % ~1e-15: and the other way round

Do not test the null spaces by norm(null(A) - null®) — that comes out around 2, not 0. Both are orthonormal bases of the same subspace, but rotated relative to each other, so comparing them entry by entry compares the bases and not the spaces. Test the subspace instead, as above.

The complete solution of \(Ax = b\)

Check solvability first, then assemble:

rank(A) == rank([A b])   % 1 means b is in C(A)
5*b(1) - 2*b(2) + b(3)   % the solvability condition from the notes: 0
xp = [-2; 0; 1; 0];      % particular solution
A*xp - b                 % zero
N  = [-3 1; 1 0; 0 -1; 0 1];   % the two special solutions, as columns
x  = xp + N*[3; -1]      % any coefficients give another solution
A*x - b                  % still zero

Backslash also works, but read the warning it prints:

x = A\b

For a rank-deficient system Matlab warns ("Rank deficient, rank = 2") and returns a basic solution, here \((0,-2/3,1,0)\) — with as many zeros as it can manage, which is neither the \(x_p\) of the notes nor the shortest solution. pinv(A)*b returns the shortest one, \((-0.143,-0.429,0.429,0.571)\).

Octave differs on both counts: it prints no warning and returns the shortest solution, i.e. the same answer as pinv. So on Octave the paragraph above describes pinv, not backslash. Either way, neither is "the" answer — the complete solution is the whole family assembled above.

The four subspaces and their right angles

r = rank(A); [m,n] = size(A);
[r, r, n-r, m-r]         % dimensions of C(A^T), C(A), N(A), N(A^T): 2 2 2 1
null(A')                 % the left null space, 3 x 1

The notes say this is spanned by \((5,-2,1)\). null normalises, so compare directions:

y = null(A'); y/y(3)     % (5, -2, 1) up to round off
y' * A                   % zero: y is perpendicular to every column

Now verify the two perpendicularity claims:

norm( orth(A')' * null(A) )    % row space  vs  null space:  0
norm( orth(A)'  * null(A') )   % column space vs left null space: 0

One-sided inverses

F = [1 2 3; 4 5 7];          % 2 x 3, full row rank
C = F'*inv(F*F');            % right inverse
F*C                          % I(2)
C*F                          % NOT I(3)

G = F';                      % 3 x 2, full column rank
B = inv(G'*G)*G';            % left inverse
B*G                          % I(2)
G*B                          % NOT I(3)

Taking \(G = F^T\) is deliberate: full row rank for \(F\) is the same statement as full column rank for \(F^T\). That is also why C*F and G*B come out equal — both are the projection onto the same two-dimensional space, which is an idea from the chapter on orthogonality arriving early.

pinv gives the same two matrices, and is what you should actually use — it does not form the inverse explicitly:

norm(pinv(F) - C), norm(pinv(G) - B)

Transformations as matrices

Differentiation and integration on polynomials, with the coefficient ordering of the notes (constant term first). Everything here is a coordinate column, never a polynomial: p below is \([p]_B\), five long because the basis \(B\) of \(P_4\) has five vectors.

Adiff = [0 1 0 0 0; 0 0 2 0 0; 0 0 0 3 0; 0 0 0 0 4];    % P4 -> P3, 4 x 5
Aint  = [0 0 0 0; 1 0 0 0; 0 1/2 0 0; 0 0 1/3 0; 0 0 0 1/4];  % P3 -> P4, 5 x 4
Adiff*Aint               % I(4): differentiate after integrating
Aint*Adiff               % NOT I(5): the constant term is lost
p = [0; 1; 2; -4; 1];    % the COORDINATES of t + 2t^2 - 4t^3 + t^4 in B
Adiff*p                  % (1, 4, -12, 4): the coordinates of the derivative, in C

Cross-check against the built-ins, which order coefficients the other way round (highest power first), so the vector must be flipped:

flip(polyder(flip(p')))  % same answer

Seeing a transformation act

The most instructive thing here is to watch a whole shape move, rather than a single vector.

S   = [0 1 1 0 0; 0 0 1 1 0];        % the unit square, as 5 points
th  = pi/6;
Rot = [cos(th) -sin(th); sin(th) cos(th)];
P   = [cos(th)^2 sin(th)*cos(th); sin(th)*cos(th) sin(th)^2];
H   = 2*P - eye(2);

plot(S(1,:), S(2,:), 'k-', 'LineWidth', 1.5); hold on; axis equal; grid on
T = Rot*S;  plot(T(1,:), T(2,:), 'b-', 'LineWidth', 1.5);   % rotated
T = P*S;    plot(T(1,:), T(2,:), 'r-', 'LineWidth', 2);     % flattened onto a line
T = H*S;    plot(T(1,:), T(2,:), 'g-', 'LineWidth', 1.5);   % reflected
legend('original','rotated','projected','reflected')

The projected square collapses to a segment — that is a null space you can see. Then check the algebra from the notes:

norm(P*P - P)            % P^2 = P
norm(H*H - eye(2))       % H^2 = I
rank(P)                  % 1
null(P)                  % the direction perpendicular to the projection line

Finally, the change of basis. Put the projection line and its perpendicular into the columns of Q and look at the same transformations in that basis:

u = [cos(th); sin(th)]; w = [-sin(th); cos(th)];
Q = [u w];
Q\P*Q                    % diag(1,0)
Q\H*Q                    % diag(1,-1)

Two awkward matrices become diagonal, only by choosing sensible axes. This is the idea behind eigenvectors, which is where the course goes next.

Home