Matrix projection blues

savecancel

Q1. In his lecture, Prof. Strang begins with the simplest case: the orthogonal projection of a vector $\mathbf{b}$ onto a one dimensional subspace $S=\mathrm{span}\langle\mathbf{a}\rangle$, $\mathbf{a}$ being a non-null vector in $\mathbb{R}^n$. The projection vector $\mathbf{p}$ should be a multiple of $\mathbf{a}$, so Prof. Strang writes it as $\mathbf{p}=x\mathbf{a}$. Here $x$ is a scalar number, not a matrix, and $\mathbf{a}$ is a vector, identified with a $n\times 1$ matrix. He wants to deduce that $\mathbf{p}=P\mathbf{b}$, where $P$ is the projection matrix

$$P=\frac{1}{\mathbf{a}^T\mathbf{a}}\mathbf{a}\mathbf{a}^T,$$ which is an $n\times n$ matrix. To this end, Prof. Strang writes the projection vector as $\mathbf{p}=\mathbf{a}\,x$. Now, you should see $\mathbf{p}$ and $\mathbf{a}$ as $n\times1$ matrices and $x$ as a one dimensional vector or a $1\times1$ matrix, so that $\mathbf{a}$ and $x$ can be multiplied. Thus, we have $$x\,\mathbf{a}=\mathbf{a}\,x.$$ This identity does not mean that $\mathbf{a}$ and $x$ are two matrices that commute. Just read the left and right sides as stated above. It is true that the notation may be a bit misleading. Perhaps it would be better to write $(x)$ instead of $x$ if $x$ should be seen as a $1\times 1$ matrix (parentheses are used to delimit matrices). For example, $$4\begin{pmatrix} 1 \\ 2 \\ 3\end{pmatrix} =\begin{pmatrix} 4 \\ 8 \\ 12\end{pmatrix} =\begin{pmatrix} 1 \\ 2 \\ 3\end{pmatrix}(4).$$

The advantage of expressing $\mathbf{p}$ as $\mathbf{a}x$ is that generalization is then possible. If $S$ is spanned by the linearly independent vectors $\mathbf{a}_1,\ldots,\mathbf{a}_k$, the orthogonal projection $\mathbf{p}$ of a vector $\mathbf{b}$ onto $S$ should be a linear combination of $\mathbf{a}_1,\ldots,\mathbf{a}_k$, that is,

$$\mathbf{p}=x_1\mathbf{a}_1+\cdots+x_k\mathbf{a}_k=AX$$ with $$A=\left(\mathbf{a}_1 \vert\ldots \vert\mathbf{a}_k\right)$$ and $$X=\begin{pmatrix}x_1 \\ \vdots \\ x_k\end{pmatrix}$$ The reasoning in Prof. Strang’s lecture would then show that $\mathbf{p}=P\mathbf{b}$, where $P$ is now the projection matrix $$P=A(A^TA)^{-1}A^T.$$

Returning to your question, you write P=X*A where, in fact, you should write $\mathbf{p}=x\mathbf{a}$ and consider $x$ being an scalar, not a matrix. Consequently, since $x$ is a scalar, from $\mathbf{a}^T(\mathbf{b}-\mathbf{p})=0$, one has

$$\mathbf{a}^T\mathbf{b}=\mathbf{a}^T\mathbf{p}=\mathbf{a}^T(x\mathbf{a})=x\mathbf{a}^T\mathbf{a}.$$ Hence $$x=\frac{\mathbf{a}^T\mathbf{b}}{\mathbf{a}^T\mathbf{a}}.$$

Q2. In your code, A.transpose()*A is a $1\times 1$ matrix. You cannot divide by a matrix. Thus you need to extract the unique element of this matrix, either with the .det() method or simply writing (A.transpose()*A)[0,0]. By the way, note that, in the video, the projection matrix is denoted by $P$, not $X$.

In SageMath, I think that it is better to use vector instead of matrix to represent $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{p}$, as shown here.