Principal component analysis can use matrix decomposition, but that is just a tool to get there.

How would you find the principal components without the use of matrix algebra?

What is the objective function (goal), and what are the constraints?

Without trying to give a full primer on PCA, from an optimization standpoint, the primary objective function is the Rayleigh quotient. The matrix that figures in the quotient is (some multiple of) the sample covariance matrix
$$\newcommand{\m}[1]{\mathbf{#1}}\newcommand{\x}{\m{x}}\newcommand{\S}{\m{S}}\newcommand{\u}{\m{u}}\newcommand{\reals}{\mathbb{R}}\newcommand{\Q}{\m{Q}}\newcommand{\L}{\boldsymbol{\Lambda}}

\S = \frac{1}{n} \sum_{i=1}^n \x_i \x_i^T = \m{X}^T \m{X} / n
$$
where each $\x_i$ is a vector of $p$ features and $\m{X}$ is the matrix such that the $i$th row is $\x_i^T$.

PCA seeks to solve a sequence of optimization problems. The first in the sequence is the unconstrained problem
$$
\begin{array}{ll}
\text{maximize} & \frac{\u^T \S \u}{\u^T\u} \;, \u \in \reals^p \> .
\end{array}
$$

Since $\u^T \u = \|\u\|_2^2 = \|\u\| \|\u\|$, the above unconstrained problem is equivalent to the constrained problem
$$
\begin{array}{ll}
\text{maximize} & \u^T \S \u \\
\text{subject to} & \u^T \u = 1 \>.
\end{array}
$$

Here is where the matrix algebra comes in. Since $\S$ is a symmetric positive semidefinite matrix (by construction!) it has an eigenvalue decomposition of the form
$$
\S = \Q \L \Q^T \>,
$$
where $\Q$ is an orthogonal matrix (so $\Q \Q^T = \m{I}$) and $\L$ is a diagonal matrix with nonnegative entries $\lambda_i$ such that $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p \geq 0$.

Hence, $\u^T \S \u = \u^T \Q \L \Q^T \u = \m{w}^T \L \m{w} = \sum_{i=1}^p \lambda_i w_i^2$. Since $\u$ is constrained in the problem to have a norm of one, then so is $\m{w}$ since $\|\m{w}\|_2 = \|\Q^T \u\|_2 = \|\u\|_2 = 1$, by virtue of $\Q$ being orthogonal. 

But, if we want to maximize the quantity $\sum_{i=1}^p \lambda_i w_i^2$ under the constraints that $\sum_{i=1}^p w_i^2 = 1$, then the best we can do is to set $\m{w} = \m{e}_1$, that is, $w_1 = 1$ and $w_i = 0$ for $i > 1$.

Now, backing out the corresponding $\u$, which is what we sought in the first place, we get that
$$
\u^\star = \Q \m{e}_1 = \m{q}_1 
$$
where $\m{q}_1$ denotes the first column of $\Q$, i.e., the eigenvector corresponding to the largest eigenvalue of $\S$. The value of the objective function is then also easily seen to be $\lambda_1$.



The remaining principal component vectors are then found by solving the sequence (indexed by $i$) of optimization problems
$$
\begin{array}{ll}
\text{maximize} & \u_i^T \S \u_i \\
\text{subject to} & \u_i^T \u_i = 1 \\
                & \u_i^T \u_j = 0 \quad \forall 1 \leq j < i\>.
\end{array}
$$
So, the problem is the same, except that we add the additional constraint that the solution must be orthogonal to all of the previous solutions in the sequence. It is not difficult to extend the argument above inductively to show that the solution of the $i$th problem is, indeed, $\m{q}_i$, the $i$th eigenvector of $\S$.



The PCA solution is also often expressed in terms of the singular value decomposition of $\m{X}$. To see why, let $\m{X} = \m{U} \m{D} \m{V}^T$. Then $n \S = \m{X}^T \m{X} = \m{V} \m{D}^2 \m{V}^T$ and so $\m{V} = \m{Q}$ (strictly speaking, up to sign flips) and $\L = \m{D}^2 / n$.

The principal components are found by projecting $\m{X}$ onto the principal component vectors. From the SVD formulation just given, it's easy to see that
$$
\m{X} \m{Q} = \m{X} \m{V} = \m{U} \m{D} \m{V}^T \m{V} = \m{U} \m{D} \> .
$$

The simplicity of representation of both the principal component vectors and the principal components themselves in terms of the SVD of the matrix of features is one reason the SVD features so prominently in some treatments of PCA.
