Coding recap

Section 6.6 Coding recap

Subsection Linked notebook

sebroc.github.io/MATH345-Course-Materials/labs/#lab-u6

. The lab checks one- and several-constraint Lagrange conditions, connects quadratic forms with maximum stretch, computes PCA from covariance and from the SVD, compares PCA with linear regression, and reads rank and the four fundamental subspaces from an SVD.

🔗

For a quick reference on arrays and shapes, matrix products, numerical checks, eigenvalue and SVD commands, and plotting, see the programming appendix sections B.2, B.4, B.9, B.10, and B.14.

🔗

Subsection Review activities

These activities interpret Python code and output for PCA and SVD: centering data, covariance matrices, principal directions, reconstruction, reduced-SVD shapes, and singular-vector identities. Reopen the definition “First principal direction” or the definition “Singular value decomposition” as needed.

🔗

Activity 6.6.1. Reading a PCA computation (U6-LO5).

The following code repeats part of the activity “Principal directions of a small data set”.

🔗

import numpy as np

np.set_printoptions(precision=3, suppress=True)

X = np.array([
    [4.0, 3.0],
    [3.0, 4.0],
    [0.0, 1.0],
    [1.0, 0.0],
])

xbar = X.mean(axis=0)
Z = X - xbar
C = Z.T @ Z / len(X)

eigvals, V = np.linalg.eigh(C)
order = np.argsort(eigvals)[::-1]
eigvals = eigvals[order]
V = V[:, order]

v1 = V[:, 0]
scores1 = Z @ v1
Zhat1 = np.outer(scores1, v1)
R1 = Z - Zhat1

(
    xbar,
    C,
    eigvals,
    V,
    Zhat1,
    np.allclose(R1 @ v1, np.zeros(len(X))),
    V.T @ C @ V,
)

Output:

🔗

(array([2., 2.]),
 array([[2.5, 2. ],
        [2. , 2.5]]),
 array([4.5, 0.5]),
 array([[ 0.707, -0.707],
        [ 0.707,  0.707]]),
 array([[ 1.5,  1.5],
        [ 1.5,  1.5],
        [-1.5, -1.5],
        [-1.5, -1.5]]),
 True,
 array([[4.5, 0. ],
        [0. , 0.5]]))

Give the shapes of X, xbar, Z, C, V, scores1, and Zhat1.
🔗

🔗
What mathematical objects are stored in xbar and C?
🔗

🔗
Why does the code reverse the order returned by np.linalg.eigh?
🔗

🔗
What do the columns of V represent?
🔗

🔗
Interpret scores1, Zhat1, and R1.
🔗

🔗
What geometric condition is checked by
🔗
```
np.allclose(R1 @ v1, np.zeros(len(X)))
```
?
🔗

🔗
Another valid numerical output could replace a column of V by its negative. What would happen to the corresponding scores and reconstructions?
🔗

🔗
Interpret the final diagonal matrix. What fraction of the total centered variation is captured by the first principal direction?
🔗

🔗

Solution.

The shapes are

\begin{equation*} X\in\R^{4\times2}, \qquad \bar{\mathbf{x}}\in\R^2, \qquad Z\in\R^{4\times2}, \end{equation*}

\begin{equation*} C\in\R^{2\times2}, \qquad V\in\R^{2\times2}, \qquad \mathbf{t}_1\in\R^4, \qquad \widehat Z_1\in\R^{4\times2}. \end{equation*}

🔗

The array xbar stores the sample mean

\begin{equation*} \bar{\mathbf{x}} = \begin{bmatrix} 2\\ 2 \end{bmatrix}, \end{equation*}

while C stores the covariance matrix

\begin{equation*} C= \begin{bmatrix} 5/2 \amp 2\\ 2 \amp 5/2 \end{bmatrix}. \end{equation*}

🔗

The command np.linalg.eigh returns the eigenvalues in increasing order. The reversal places the largest eigenvalue first. The first column of \(V\) is therefore a first principal direction, and the second column is a second principal direction.

🔗

The vector scores1 stores

\begin{equation*} \mathbf{t}_1=Z\mathbf{v}_1. \end{equation*}

The rows of Zhat1 are the projections of the centered data vectors onto

\begin{equation*} \spans\{\mathbf{v}_1\}. \end{equation*}

The matrix R1 stores the corresponding reconstruction residuals.

🔗

The product R1 @ v1 computes the dot product of each residual row with \(\mathbf{v}_1\text{.}\) The value True checks that every reconstruction residual is numerically perpendicular to the first principal direction.

🔗

Replacing \(\mathbf{v}_1\) by \(-\mathbf{v}_1\) reverses the signs of its scores. However,

\begin{equation*} (-t_i)(-\mathbf{v}_1)=t_i\mathbf{v}_1, \end{equation*}

so the projected points and reconstruction residuals do not change. The principal axis is the same.

🔗

The final output says that covariance in the principal-coordinate system is

\begin{equation*} V^TCV = \begin{bmatrix} 4.5 \amp 0\\ 0 \amp 0.5 \end{bmatrix}. \end{equation*}

The diagonal entries are the captured variations in the two principal directions. The first direction captures

\begin{equation*} \frac{4.5}{4.5+0.5} = \frac9{10} \end{equation*}

of the total centered variation.

🔗

Activity 6.6.2. Reading SVD shapes (U6-LO6).

The following code computes a reduced SVD.

🔗

import numpy as np

A = np.array([[ 1., -1.],
              [-2.,  2.],
              [ 2., -2.]])

U, s, Vt = np.linalg.svd(A, full_matrices=False)

A.shape, U.shape, s.shape, Vt.shape

Output:

🔗

((3, 2), (3, 2), (2,), (2, 2))

What is the shape of \(A\text{?}\)
🔗

🔗
How many singular values are returned?
🔗

🔗
Why is s.shape equal to (2,) rather than (2, 2)?
🔗

🔗
What matrix does Vt represent?
🔗

🔗
What does full_matrices=False do in this example?
🔗

🔗

Solution.

The matrix \(A\) is \(3\times 2\text{.}\) The reduced SVD returns two singular values because

\begin{equation*} \min(3,2)=2. \end{equation*}

The array s stores only the diagonal entries of \(\Sigma\text{,}\) not the full diagonal matrix, so its shape is (2,). The array Vt represents \(V^T\text{.}\) With full_matrices=False, NumPy returns the reduced SVD shapes needed for reconstruction:

\begin{equation*} U\in\mathbb R^{3\times 2}, \qquad s\in\mathbb R^2, \qquad V^T\in\mathbb R^{2\times 2}. \end{equation*}

🔗

Activity 6.6.3. Checking a singular-vector identity in code (U6-LO6).

For an SVD

\begin{equation*} A=U\Sigma V^T, \end{equation*}

the identity

\begin{equation*} A\mathbf{v}_i=\sigma_i\mathbf{u}_i \end{equation*}

is one of the main ways to read the factors. The following code checks this identity for a diagonal stretch.

🔗

A = np.array([[3., 0.],
              [0., 1.]])

U, s, Vt = np.linalg.svd(A, full_matrices=False)

i = 0
v = Vt.T[:, i]
lhs = A @ v
rhs = s[i] * U[:, i]

lhs, rhs, np.allclose(lhs, rhs)

Output:

🔗

(array([3., 0.]), array([3., 0.]), True)

Which vector is stored in v?
🔗

🔗
What does lhs compute?
🔗

🔗
What does rhs compute?
🔗

🔗
Why does np.allclose(lhs, rhs) return True?
🔗

🔗
Which stretch factor appears in this computation?
🔗

🔗

Solution.

The vector v is the first right singular vector \(\mathbf{v}_1\text{.}\) The array lhs computes \(A\mathbf{v}_1\text{.}\) The array rhs computes \(\sigma_1\mathbf{u}_1\text{.}\) The value True means that the numerical computation agrees with

\begin{equation*} A\mathbf{v}_1=\sigma_1\mathbf{u}_1. \end{equation*}

The stretch factor is

\begin{equation*} \sigma_1=3. \end{equation*}

🔗

Prev Top Next