Linear Algebra Review and SciPy Basics

Definitions and notation in this review follow Section 2.2 Vectors and Matrices in Biegler (2010), with the invertibility discussion cross-checked against Lay, Linear Algebra and Its Applications, 3rd ed.

Recommended Links

• NumPy Documentation - Main Page

• SciPy Documentation - Main Page

• SciPy Linear Algebra - Tutorial

• SciPy Linear Algebra - API

• SciPy Lecture Notes (Especially 1.1 - 1.6 for everyone and 2.1 - 2.2 for advanced users)

Tip: We will mostly use SciPy for linear algebra. (It has more sophisticated capabilities than NumPy.) We will use NumPy if/when SciPy does not offer a specific command.

This guide supplies definitions, proof details, and extended calculations for the lecture. Use the sections on spaces and rank, worked factorizations, SVD, and conditioning alongside the shorter handout. Numerical Errors computes roundoff and error measurements. The scaling example returns scaled variables to their original units.

Notation and linear maps¶

Scalars and scalar components are plain, vectors are bold lower case, and matrices are bold upper case: $\alpha$, $x_i$, $\mathbf{x}$, $A_{ij}$, and $\mathbf{A}$. The matrix $\mathbf{A}\in\mathbb{R}^{n\times m}$ represents the linear map $T:\mathbb{R}^m\to\mathbb{R}^n$, $T(\mathbf{x})=\mathbf{Ax}$. Linearity means

for every $\mathbf{x},\mathbf{y}\in\mathbb{R}^m$ and $\alpha,\beta\in\mathbb{R}$. The columns of $\mathbf{A}$ are the images of the coordinate vectors: $\mathbf{A}\mathbf{e}_j=\mathbf{a}_j$.

Matrix operations¶

For $\mathbf{B}\in\mathbb{R}^{m\times p}$, $\mathbf{C}=\mathbf{AB}\in\mathbb{R}^{n\times p}$ has entries

The inner dimensions must match. Multiplication composes maps: first apply $\mathbf{B}$, then $\mathbf{A}$. If $\mathbf{A}$ is $4\times3$ and $\mathbf{B}$ is $3\times7$, $\mathbf{AB}$ is $4\times7$; $\mathbf{BA}$ is undefined because its inner dimensions are 7 and 4.

The transpose interchanges rows and columns, $(\mathbf{A}^\mathsf{T})_{ij}=A_{ji}$. A square matrix is symmetric when $\mathbf{A}=\mathbf{A}^\mathsf{T}$ and diagonal when $A_{ij}=0$ for $i\ne j$. The identity has $I_{ii}=1$ and all off-diagonal entries zero.

For the illustrated $\mathbf{A}=\left[\begin{smallmatrix}1.5&.5\\.5&1\end{smallmatrix}\right]$, the area changes by $\det\mathbf{A}=1.25$. A negative determinant reverses orientation; a zero determinant collapses area to zero.

Video: 3Blue1Brown, Linear transformations and matrices and Matrix multiplication as composition.

Activity

• Calculate $\mathbf{A}_o^{T}$

• Create a square matrix and check if it is symmetric

• Create a 3x3 identity matrix

Determinant¶

For a square matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$, expansion along any row $i$ gives

where $A_{ij}$ is the entry in row $i$ and column $j$, and $\bar A_{ij}$ is the determinant of the $(n-1)\times(n-1)$ matrix obtained by deleting row $i$ and column $j$ from $\mathbf{A}$. Useful properties are

The exponent $n$ in the scaling law matters: multiplying the whole matrix by $\alpha$ scales each of its $n$ rows. A square matrix is singular exactly when $\det(\mathbf{A})=0$.

Activity

• Verify the properties listed under Some properties for the determinant include

• Are $\mathbf{A}_d$, $\mathbf{B}_d$ or $\mathbf{C}_d$ singular?

Another Example (for at home)

Expand along the first row of

Then

Span, independence, and rank¶

A linear combination of $\mathbf{a}_1,\ldots,\mathbf{a}_m\in\mathbb{R}^n$ is $\sum_{j=1}^m\alpha_j\mathbf{a}_j$, with real coefficients. Their span is the set of all such combinations. They are linearly independent if $\sum_j\alpha_j\mathbf{a}_j=\mathbf{0}$ implies $\alpha_j=0$ for every $j$; otherwise they are dependent.

A subspace $S\subseteq\mathbb{R}^n$ contains $\mathbf{0}$ and is closed under linear combinations: $\alpha\mathbf{x}+\beta\mathbf{y}\in S$ for all $\mathbf{x},\mathbf{y}\in S$ and all real $\alpha,\beta$. A basis of $S$ is a linearly independent spanning set. Its number of vectors is $\dim S$. The zero subspace has the empty basis and dimension zero.

For $\mathbf{A}\in\mathbb{R}^{n\times m}$, the range or column space is

The null space is $\operatorname{Null}(\mathbf{A})=\{\mathbf{x}\in\mathbb{R}^m:\mathbf{Ax}=\mathbf{0}\}$. Both are subspaces, but they belong to different spaces: range in $\mathbb{R}^n$, null space in $\mathbb{R}^m$.

The rank $r(\mathbf{A})$ is the dimension of the column space, equivalently the row space or the size of the largest square submatrix with nonzero determinant. Full column rank means $r=m$; full row rank means $r=n$. Full rank means $r=\min(n,m)$; otherwise the matrix is rank deficient.

Rank-nullity: complete a basis¶

Let $\mathbf{v}_1,\ldots,\mathbf{v}_k$ be a basis of the null space, and extend it to a basis $\mathbf{v}_1,\ldots,\mathbf{v}_m$ of the domain. Every $\mathbf{Ax}$ is a linear combination of $\mathbf{A}\mathbf{v}_{k+1},\ldots,\mathbf{A}\mathbf{v}_m$, because the first $k$ vectors map to zero. These remaining images are independent: any combination mapping to zero would be both a combination of the last basis vectors and a member of the span of the first $k$. Independence of the whole basis forces every coefficient to vanish. Thus

For the illustrated matrix $\left[\begin{smallmatrix}.5&1\\.5&1\end{smallmatrix}\right]$, the range is $\operatorname{span}\{(1,1)^\mathsf{T}\}$ and the null space is $\operatorname{span}\{(2,-1)^\mathsf{T}\}$. Both have dimension one. If $\mathbf{x}_0$ solves $\mathbf{Ax}=\mathbf{b}$, every solution is $\mathbf{x}_0+\mathbf{z}$ with $\mathbf{z}\in\operatorname{Null}(\mathbf{A})$. If $\mathbf{b}$ is outside the range, there is no solution.

Video: 3Blue1Brown, Linear combinations, span, and basis vectors, and Inverse matrices, column space, and null space.

Activity. Predict the ranks of $\mathbf{A}_o^\mathsf{T}$, $\mathbf{B}_o^\mathsf{T}$, $\mathbf{A}_d$, $\mathbf{B}_d$, and $\mathbf{C}_d$. Compare with np.linalg.matrix_rank and explain any rank deficiency in terms of dependent columns.

Inverse¶

A square matrix $\mathbf{A}$ is invertible when there is a matrix $\mathbf{A}^{-1}$ satisfying $\mathbf{AA}^{-1}=\mathbf{A}^{-1}\mathbf{A}=\mathbf{I}$. Then $\mathbf{Ax}=\mathbf{b}$ has the unique solution $\mathbf{x}=\mathbf{A}^{-1}\mathbf{b}$, although numerical software should solve the system directly rather than form the inverse. Also, $\det(\mathbf{A}^{-1})=1/\det(\mathbf{A})$. An orthogonal matrix $\mathbf{Q}$ satisfies $\mathbf{Q}^{-1}=\mathbf{Q}^\mathsf{T}$.

Activity

1. Calculate the inverse of the square matrix $\mathbf{A}_d$ and verify that $\mathbf{A}_d^{-1} \mathbf{A}_d = \mathbf{I}$. (The rectangular matrix $\mathbf{A}_o$ has no two-sided inverse.)

2. Verify that $\det(\mathbf{A}_d^{-1}) = 1/\det(\mathbf{A}_d)$

3. Is $\mathbf{Q}$ an orthogonal matrix?

If $\det\mathbf{A}=0$ and an inverse existed, multiplicativity of the determinant would give $1=\det\mathbf{I}=\det(\mathbf{AA}^{-1})=\det\mathbf{A}\det\mathbf{A}^{-1}=0$, a contradiction. The Invertible Matrix Theorem collects the converse and equivalent conditions.

For dense unstructured matrices, factorizing and solving usually takes less work than explicitly forming an inverse and multiplying. Leading operation counts for a single right-hand side are approximately $2n^3/3$ for elimination versus $2n^3$ for an inverse. Accuracy also depends on the algorithm and data; avoid treating a computed inverse as a universal cure for a difficult solve. See Higham (2002), §14.1, p.260.

A sparse matrix can have a dense inverse, increasing storage as well as arithmetic. For an explicit example, let $\mathbf N$ have ones immediately above the diagonal and zeros elsewhere. Then $\mathbf N^n=\mathbf0$ and

Multiplication telescopes to $\mathbf I$. The original matrix has $2n-1$ nonzeros, while its inverse has $n(n+1)/2$ nonzeros. This example explains the storage risk; sparse factorization can also create fill and needs an appropriate ordering.

Solving Linear Systems¶

Consider the linear system $\mathbf{A}_l \mathbf{x} = b_l$

Explicit Inverse¶

Calculate $\mathbf{x}$ by explicitly using $\mathbf{A}_l^{-1}$. Hint: Use linalg.inv().

LU Decomposition¶

Perform LU decomposition on $\mathbf{A}_l$. What structures do the $\mathbf{P}$, $\mathbf{L}$, and $\mathbf{U}$ matrices have?

SciPy returns $\mathbf{A}=\mathbf{PLU}$, so the solve is $\mathbf{Ly}=\mathbf{P}^\mathsf{T}\mathbf{b}$ followed by $\mathbf{Ux}=\mathbf{y}$. A permutation is orthogonal: $\mathbf{P}^{-1}=\mathbf{P}^\mathsf{T}$. It need not be symmetric. Partial pivoting interchanges rows to use a largest-magnitude entry in the active column as the pivot. This avoids division by a smaller available pivot; it does not guarantee small element growth for every matrix.

Reuse the factors¶

For a dense $n\times n$ system, factorization costs approximately $2n^3/3$ operations, and the two triangular solves cost approximately $2n^2$ per right-hand side. With $q$ right-hand sides, refactorizing each time costs $q(2n^3/3+2n^2)$; reusing the factors costs $2n^3/3+2qn^2$. For $n=1000$ and $q=200$, these are $1.3373\cdot10^{11}$ and $1.0667\cdot10^9$, a ratio of about 125.4. These leading arithmetic counts omit lower-order work and do not predict an exact wall-clock speedup.

Doolittle elimination chooses unit diagonal entries in $\mathbf{L}$. See Higham (2002), §9.2, pp.162–163, and §9.4, pp.166–167 for growth and backward error. The worked examples show the handout’s small LU system without the random matrices used below.

Is P orthogonal?¶

Yes! We see that $\mathbf{P}^T = \mathbf{P}^{-1}$ and $\mathbf{P}^T \cdot \mathbf{P} = \mathbf{I}$.

MATLAB¶

Defines LU decomposition as follows:

Consider

Using MATLAB LU definition

Substitute into linear system:

Let $\mathbf{y} = \mathbf{U} \cdot \mathbf{x}$ and substitute.

Step 1. Solve $\mathbf{L} \cdot \mathbf{y} = \mathbf{P} \cdot \mathbf{b}$ for $\mathbf{y}$

Step 2. Solve $\mathbf{U} \cdot \mathbf{x} = \mathbf{y}$ for $\mathbf{x}$

The next cell deliberately applies $\mathbf{Pb}$ to the SciPy factors computed above. This mixes the two conventions and generally fails the residual check. MATLAB’s own permutation is the transpose of SciPy’s for the same factors. Use $\mathbf{P}^\mathsf{T}\mathbf{b}$ in the SciPy solve that follows.

SciPy¶

Defines LU decomposition as:

Consider

SciPy LU definition

Substitute into linear system:

Let $\mathbf{y} = \mathbf{U} \cdot \mathbf{x}$ and substitute.

Step 1. Solve $\mathbf{L} \cdot \mathbf{y} = \mathbf{P}^T \cdot \mathbf{b}$ for $\mathbf{y}$

Step 2. Solve $\mathbf{U} \cdot \mathbf{x} = \mathbf{y}$ for $\mathbf{x}$

What is the difference? A transpose.

Verify our answer with linalg.solve¶

Solve the linear system using linalg.solve

Invertible Matrix Theorem¶

For a square $n\times n$ matrix $\mathbf{A}$, the following statements are equivalent:

• $\mathbf{A}$ is invertible; its row-reduced form is $\mathbf{I}$; and it has $n$ pivots.

• $\mathbf{Ax}=\mathbf{0}$ has only the trivial solution, so $\operatorname{Nul}(\mathbf{A})=\{\mathbf{0}\}$ and $\dim\operatorname{Nul}(\mathbf{A})=0$.

• The columns of $\mathbf{A}$ are linearly independent, form a basis of $\mathbb{R}^n$, and span $\mathbb{R}^n$; equivalently $\operatorname{rank}(\mathbf{A})=n$.

• The map $\mathbf{x}\mapsto \mathbf{Ax}$ is both one-to-one and onto; therefore $\mathbf{Ax}=\mathbf{b}$ has a solution for every $\mathbf{b}\in\mathbb{R}^n$.

• A left or right inverse exists; $\mathbf{A}^\mathsf{T}$ is invertible; $\det(\mathbf{A})\ne0$; and zero is not an eigenvalue.

• The row space is $\mathbb{R}^n$, $(\operatorname{Col}\mathbf{A})^\perp=\{\mathbf{0}\}$, $(\operatorname{Nul}\mathbf{A})^\perp=\mathbb{R}^n$, and $\mathbf{A}$ has $n$ nonzero singular values.

These equivalences consolidate Lay’s Invertible Matrix Theorem into one searchable list.

Eigenvectors and Eigenvalues¶

For square $\mathbf{A}$, a nonzero vector $\mathbf{v}$ is an eigenvector with eigenvalue $\lambda$ when $\mathbf{Av}=\lambda \mathbf{v}$. Thus $\mathbf{A}-\lambda \mathbf{I}$ is singular and eigenvalues are roots of $\det(\mathbf{A}-\lambda \mathbf{I})=0$. Counting algebraic multiplicity, $\det(\mathbf{A})=\prod_i\lambda_i$.

If $\mathbf{A}$ is real and symmetric, the spectral theorem gives $\mathbf{A}=\mathbf{V}\boldsymbol{\Lambda} \mathbf{V}^\mathsf{T}$ with real eigenvalues and orthonormal eigenvectors. Only under this symmetry assumption can definiteness be read directly from eigenvalue signs: all positive means positive definite, all nonnegative means positive semidefinite, and mixed signs mean indefinite.

A nonsymmetric real matrix may have complex eigenvalues, and its eigenvalues do not in general classify the quadratic form $\mathbf{x}^\mathsf{T}\mathbf{Ax}$. That form is governed by the symmetric part $(\mathbf{A}+\mathbf{A}^\mathsf{T})/2$. For rectangular matrices, singular values $\sigma_i\ge0$—the square roots of the eigenvalues of $\mathbf{A}^\mathsf{T}\mathbf{A}$—play the analogous role.

A nonsymmetric worked example¶

For $\mathbf{A}=\left[\begin{smallmatrix}0&1\\-2&-3\end{smallmatrix}\right]$,

For $\lambda_1=-1$, the equations in $(\mathbf{A}+\mathbf{I})\mathbf{v}=\mathbf{0}$ are multiples of $v_1+v_2=0$. The matrix has rank one, and eigenvectors are $k(1,-1)^\mathsf{T}$ with $k\ne0$. For $\lambda_2=-2$ the eigenvectors are $k(1,-2)^\mathsf{T}$ with $k\ne0$. These eigendirections are not orthogonal; symmetry is the hypothesis that guarantees an orthonormal eigenbasis. Counting algebraic multiplicity, $\operatorname{tr}\mathbf{A}=\sum_i\lambda_i$ and $\det\mathbf{A}=\prod_i\lambda_i$. For nonsingular $\mathbf{A}$, the eigenpairs of $\mathbf{A}^{-1}$ are $(1/\lambda_i,\mathbf{v}_i)$.

Video: 3Blue1Brown, Eigenvectors and eigenvalues.

Definiteness and the spectral theorem¶

For real symmetric $\mathbf{A}$ and every nonzero real $\mathbf{x}$, positive definite means $\mathbf{x}^\mathsf{T}\mathbf{Ax}>0$; positive semidefinite means $\mathbf{x}^\mathsf{T}\mathbf{Ax}\ge0$. Reverse the signs for negative definite or negative semidefinite. Indefinite means the form takes both positive and negative values. Singular means the matrix is not invertible; it is a different property from definiteness.

Write $\mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^\mathsf{T}$ with orthogonal $\mathbf{V}$ and set $\mathbf{z}=\mathbf{V}^\mathsf{T}\mathbf{x}$. Then

This proves the eigenvalue sign tests. Choosing $\mathbf{x}=\mathbf{v}_i$ proves necessity; summing nonnegative squared coordinates proves sufficiency. A symmetric semidefinite matrix can be singular: $\left[\begin{smallmatrix}1&1\\1&1\end{smallmatrix}\right]$ has eigenpairs $(2,(1,1)^\mathsf{T})$ and $(0,(1,-1)^\mathsf{T})$. Its quadratic form is $(x_1+x_2)^2$, so it is positive semidefinite and has a nontrivial null space.

For symmetric nonsingular $\mathbf A=\mathbf V\boldsymbol\Lambda\mathbf V^{\mathsf T}$,

Multiplying the two factorizations gives $\mathbf I$. Thus inversion preserves each eigenvector and reciprocates its eigenvalue. Keep each eigenvalue paired with its direction when reading a quadratic form or its inverse.

Level sets of a quadratic form, as one picture¶

Eigenvalues and eigenvectors of a symmetric matrix have a picture: they are the principal axes of the level sets of $q(\mathbf{x}) = \tfrac12 \mathbf{x}^\mathsf{T} \mathbf{A} \mathbf{x}$. In the rotated coordinates $\mathbf{z} = \mathbf{V}^\mathsf{T} \mathbf{x}$ the level set is $\lambda_1 z_1^2 + \lambda_2 z_2^2 = 2q$, so the semi-axis along $v_i$ has length $\sqrt{2q/\lambda_i} \propto 1/\sqrt{\lambda_i}$ --- a large eigenvalue means steep, hence short.

The right-hand panel is the same picture with a condition number of 100. Here $K(\mathbf{A})=\lambda_{\max}/\lambda_{\min}$ is the condition number of a symmetric positive definite matrix (defined properly in the Condition Number section below). Note the axis ratio is $\sqrt{K(\mathbf{A})}$, not $K(\mathbf{A})$. The matrices in Four Factorizations, Four Failures below have $K(\mathbf{A})$ from about 10⁸ to 10¹⁸, and for those the ellipse is not drawable.

This figure is rendered from figures/plots/quadratic-form-level-sets.py and is the same image printed in the course pack, so what you annotate in lecture and what you see here cannot drift apart. The cells below are the live version.

Activity

• Do you expect this matrix to have negative eigenvalues? Choose: Yes / No / Cannot tell from inspection.

• Calculate the eigenvalues.

• The eigenvalues of this nonsymmetric matrix are positive. Does that prove $\mathbf{x}^\mathsf{T}\mathbf{A}_e \mathbf{x}>0$ for every nonzero $\mathbf{x}$?

• Form the symmetric part $\mathbf{H}=(\mathbf{A}_e+\mathbf{A}_e^\mathsf{T})/2$ and use its eigenvalues to classify the quadratic form.

• Calculate $\det(\mathbf{A}_e)$ using only the eigenvalues of $\mathbf{A}_e$.

Singular value decomposition¶

For $\mathbf{A}\in\mathbb{R}^{m\times n}$ with $m\ge n$, a thin SVD is

where $\mathbf{U}^\mathsf{T}\mathbf{U}=\mathbf{I}_n$, $\mathbf{V}^\mathsf{T}\mathbf{V}=\mathbf{VV}^\mathsf{T}=\mathbf{I}_n$, and $\boldsymbol{\Sigma}=\operatorname{diag}(\sigma_1,\ldots,\sigma_n)$ with $\sigma_1\ge\cdots\ge\sigma_n\ge0$. The columns of $\mathbf{U}$ and $\mathbf{V}$ are left and right singular vectors. For a rectangular thin $\mathbf{U}$, $\mathbf{UU}^\mathsf{T}$ is a projection, not the identity on $\mathbb{R}^m$. For $m<n$, apply the construction to $\mathbf{A}^\mathsf{T}$ or use a full SVD with rectangular $\boldsymbol{\Sigma}$.

Orthogonal factors preserve 2-norm lengths. Along $\mathbf{v}_i$, $\mathbf{A}\mathbf{v}_i=\sigma_i\mathbf{u}_i$; a zero singular value sends that direction to zero. The number of positive singular values is the exact rank. Numerical rank requires a scale-dependent tolerance: compare each singular value with a tolerance relative to $\sigma_1$, taking data accuracy into account.

Connection with a symmetric eigenproblem¶

Substitute the factorization and use $\mathbf{U}^\mathsf{T}\mathbf{U}=\mathbf{I}$:

Conversely, diagonalize the positive semidefinite $\mathbf{A}^\mathsf{T}\mathbf{A}$ and set $\mathbf{u}_i=\mathbf{Av}_i/\sigma_i$ for each positive $\sigma_i$. Then

Complete these left singular vectors to $n$ orthonormal columns if the rank is less than $n$; never divide by a zero singular value. The remaining right singular vectors lie in the null space. For full column rank, the compact formula is $\mathbf{U}=\mathbf{AV}\boldsymbol{\Sigma}^{-1}$.

This is an algebraic connection, not a recommendation to compute an SVD by forming $\mathbf{A}^\mathsf{T}\mathbf{A}$. For full column rank, its 2-norm condition number is $K_2(\mathbf{A})^2$, which can lose information about small singular values. Use scipy.linalg.svd directly. It returns Vh, meaning $\mathbf{V}^\mathsf{T}$ for real data; full_matrices=False requests the thin factors.

The pseudoinverse is $\mathbf{A}^{\dagger}=\mathbf{V}\boldsymbol{\Sigma}^{\dagger}\mathbf{U}^\mathsf{T}$, replacing positive singular values by their reciprocals and zero ones by zero. It is a two-sided inverse only for a square nonsingular matrix. Small positive singular values amplify perturbations in the corresponding directions.

See Nocedal & Wright (2006), Appendix A, p.604, Eq.(A.15), for SVD. Activity: compute the SVD of $\mathbf{A}_l$ and $\mathbf{C}_d$, reconstruct each matrix, and use the singular values to discuss its rank.

Vector and Matrix Norms¶

A vector norm $||\mathbf{x}||$ measures the length of a vector. It maps $\mathbb{R}^n$ into the nonnegative reals and satisfies three properties:

The last one is the triangle inequality. The $p$-norm is $||\mathbf{x}||_p = \left( \sum_{i=1}^n |x_i|^p \right)^{1/p}$ for $p \geq 1$, and three cases matter here:

• $||\mathbf{x}||_1 = \sum_{i=1}^n |x_i|$ --- the 1-norm

• $||\mathbf{x}||_2 = \left( \sum_{i=1}^n x_i^2 \right)^{1/2} = \sqrt{\mathbf{x}^T \mathbf{x}}$ --- the Euclidean norm

• $||\mathbf{x}||_{\infty} = \max_i |x_i|$ --- the max norm

Also $|\mathbf{x}^T \mathbf{y}| \leq ||\mathbf{x}||_p \cdot ||\mathbf{y}||_q$ for $p,q\in[1,\infty]$ with $1/p + 1/q = 1$ (using $1/\infty=0$); the case $p = q = 2$ is the Cauchy-Schwarz inequality.

Induced matrix norms are defined by vector norms:

The rules for the three induced norms are

• $||\mathbf{A}||_1 = \max_j \sum_i |A_{ij}|$ --- largest column sum

• $||\mathbf{A}||_{\infty} = \max_i \sum_j |A_{ij}|$ --- largest row sum

• $||A||_2 =$ the largest singular value of $\mathbf{A}$

Finally, the induced norms and the Frobenius norm $||\mathbf{A}||_F = \left( \sum_i \sum_j A_{ij}^2 \right)^{1/2}$ all satisfy

for any two compatible matrices $\mathbf{A}$ and $\mathbf{B}$. This is submultiplicativity, and it is what the condition number in the next section rests on.

Definition 2.7, pp. 21-22 in Biegler (2010).

The induced 2-norm is $\sigma_1$ because $\|\mathbf{Ax}\|_2^2=\sum_i\sigma_i^2z_i^2\le\sigma_1^2\|\mathbf{z}\|_2^2$ for $\mathbf{z}=\mathbf{V}^\mathsf{T}\mathbf{x}$; equality holds along a leading right singular vector. Submultiplicativity follows by applying the defining inequality twice: $\|\mathbf{ABx}\|\le\|\mathbf{A}\|\|\mathbf{Bx}\|\le\|\mathbf{A}\|\|\mathbf{B}\|\|\mathbf{x}\|$, then taking the maximum over unit vectors.

For continuous finite-dimensional vector- or matrix-valued $\mathbf{F}$ on $[a,b]$ with $a\le b$,

Apply the triangle inequality and homogeneity to each Riemann sum, then pass to the limit. Continuity ensures the sums and their norms converge. Combined with submultiplicativity, this bounds expressions such as $\|\mathbf{A}(\int\mathbf{F})\mathbf{v}\|$ by $\|\mathbf{A}\|\int\|\mathbf{F}\|\,\|\mathbf{v}\|$.

Activity

• Verify the vector norm properties hold for the p-2 norm using $\mathbf{x}$ and $\mathbf{y}$. Hint: linalg.norm()

• Calculate $||\mathbf{A}_l||_1$, $||\mathbf{A}_l||_{\infty}$ and $||\mathbf{A}_l||_2$ using both the rules given above and SciPy

Condition number and residual bounds¶

For nonsingular square $\mathbf{A}$ and an induced norm,

In the 2-norm, $K_2(\mathbf{A})=\sigma_{\max}/\sigma_{\min}$. For real symmetric nonsingular matrices this is $\max_i|\lambda_i|/\min_i|\lambda_i|$; eigenvalue magnitudes matter. Eigenvalues $2,-4,8$ give $K_2=4$. For rectangular full-column-rank matrices the singular-value ratio measures sensitivity of the map on its domain, but there is no square two-sided inverse.

Perturb only the right-hand side¶

Let $\mathbf{Ax}=\mathbf{b}\ne\mathbf{0}$ and $\mathbf{A}(\mathbf{x}+\Delta\mathbf{x})=\mathbf{b}+\Delta\mathbf{b}$. Subtract to get $\Delta\mathbf{x}=\mathbf{A}^{-1}\Delta\mathbf{b}$. Also $\|\mathbf{b}\|\le\|\mathbf{A}\|\|\mathbf{x}\|$. Therefore

The condition number is a worst-case amplification, not a prediction for every perturbation. With $\mathbf{A}=\operatorname{diag}(1,\varepsilon)$ and $\mathbf{b}=(1,0)^\mathsf{T}$, a perturbation $(0,\delta)^\mathsf{T}$ gives $\Delta\mathbf{x}=(0,\delta/\varepsilon)^\mathsf{T}$ and attains $K_2=1/\varepsilon$ for $0<\varepsilon\le1$. Perturbing along the first coordinate is amplified only by one.

For an approximate solution $\hat{\mathbf{x}}$, define the exact residual of the stored data $\mathbf{r}=\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}$. Then $\mathbf{A}(\mathbf{x}-\hat{\mathbf{x}})=\mathbf{r}$ and the same argument gives

A relative residual 10^-12 and condition number $4.5\cdot10^{14}$ give the bound 450. This bound does not certify an accurate solution; it also does not establish that the actual forward error is large. Computing the residual in floating point introduces another error. The Numerical Errors experiment evaluates the residual and reference solution in high precision and separates rounding of the input data from error in the solve.

Backward error and stability¶

In the 2-norm, allowing perturbations only in $\mathbf{A}$ gives

for nonzero $\hat{\mathbf{x}}$. The lower bound follows from $\mathbf{r}=\Delta\mathbf{A}\hat{\mathbf{x}}$; equality is attained by $\Delta\mathbf{A}=\mathbf{r}\hat{\mathbf{x}}^\mathsf{T}/\|\hat{\mathbf{x}}\|_2^2$. If $K_2\eta_A<1$, the finite perturbation bound is

Indeed, $\|\mathbf{x}-\hat{\mathbf{x}}\|\le K_2\eta_A\|\hat{\mathbf{x}}\|\le K_2\eta_A(\|\mathbf{x}\|+\|\mathbf{x}-\hat{\mathbf{x}}\|)$; rearrange. The simpler $K_2\eta_A$ expression is a first-order rule when that product is small. If perturbations are allowed in both $\mathbf{A}$ and $\mathbf{b}$, the normwise backward error instead has denominator $\|\mathbf{A}\|\|\hat{\mathbf{x}}\|+\|\mathbf{b}\|$. State which model is used. See Higham (2002), Theorems 7.1–7.2, pp.120–121.

Conditioning describes sensitivity of a problem. Backward stability describes an algorithm’s error over its intended class of inputs. A small measured backward error on one solve is evidence about that solve; it does not prove that the algorithm is backward stable for all inputs.

Activity. Estimate $K_2(\mathbf{A}_l)$ from its SVD and compare with np.linalg.cond. Choose the perturbation direction that produces the greatest change in the answer and explain why.

Four factorizations, four numerical diagnostics¶

These examples isolate element growth in LU, loss of positive definiteness near the precision limit, loss of orthogonality in QR, and the difficulty of computing inertia near zero eigenvalues. Each table is computed by its cell. Exact digits and the first reported failure can depend on NumPy, SciPy, and the BLAS/LAPACK build. Read the actual output of your run.

The Cholesky cell catches an intended LinAlgError if the factorization breaks down. The inertia cell compares floating-point reports with an exact-rational calculation for the stored entries. A returned factorization still needs its defining properties checked: a small reconstruction residual alone does not certify orthogonality, definiteness, or correct inertia.

D1. LU: element growth on a moderately conditioned problem¶

Wilkinson’s matrix has 1 on the diagonal, -1 below it, and a final column of ones. With the pivot tie choices used here, partial pivoting makes no interchanges, and entries grow as $2^{n-1}$ during elimination. The table checks the permutation and growth against this exact-arithmetic pattern. Its reported 2-norm condition numbers stay moderate, while the error can become substantial.

The reference vector x_true is chosen first, and b = A @ x_true is rounded. Consequently the reported error includes that assembly error as well as error in the solve. The residual column is normalized by the Frobenius norm of the matrix and the 2-norm of the computed solution, as shown in the code. It is not $\|\mathbf{r}\|_2/\|\mathbf{b}\|_2$.

For this run, inspect P == I ? and compare measured growth with $2^{n-1}$. Ties are relevant: being a largest-magnitude candidate does not force every pivoting implementation to make the same choice. Large element growth weakens the usual backward-error guarantee for elimination; see Higham (2002), §9.4, pp.166–167.

Activity. Compare the error, growth, and condition number columns. Which column identifies a sensitive problem, and which identifies the risk of rounding amplification during elimination? Change the random seed and identify which quantities change.

To use a forward-error bound with this table, first match both the reference problem and the normalization. Its printed residual measure $\rho_F$ converts to the matrix-only 2-norm backward error by $\eta_A=\rho_F\|\mathbf{A}\|_F/\|\mathbf{A}\|_2$, ignoring error in evaluating the residual. The finite bound $K_2\eta_A/(1-K_2\eta_A)$ requires $K_2\eta_A<1$ and compares with the exact solution of the stored system, not automatically with x_true. The high-precision experiment on Numerical Errors makes those distinctions explicitly.

D2. Cholesky near the precision limit¶

The exact Hilbert matrix $H_{ij}=1/(i+j-1)$ is symmetric positive definite for every $n$. For any nonzero $\mathbf{z}$,

because a nonzero polynomial cannot vanish throughout the interval. A matrix stored in binary floating point has rounded entries and need not retain this property when the smallest eigenvalues are comparable to roundoff. Computing eigenvalues or a Cholesky factor introduces further rounding. The cell catches breakdown and prints it as a result.

Read the smallest computed eigenvalue and Cholesky status on each row. A negative computed eigenvalue or a Cholesky breakdown does not refute the proof for the exact rational Hilbert matrix. It concerns the stored matrix and floating-point calculation. Conversely, a successful factorization is not a high-accuracy certificate for the exact mathematical data. The numerical condition-number estimate can itself be unreliable near the precision limit.

Activity. A solver uses Cholesky breakdown to trigger a Hessian correction. State what this test reports in finite precision and what extra claim would require error analysis. Extend the loop and compare the first failure on your machine with the displayed output.

D3. QR: reconstruction and orthogonality are separate checks¶

The three columns below differ only in entries of size 10^-8. Classical Gram–Schmidt, modified Gram–Schmidt, and Householder QR can all give small reconstruction errors $\|\mathbf{QR}-\mathbf{A}\|_2$. A QR factorization also requires $\mathbf{Q}^\mathsf{T}\mathbf{Q}=\mathbf{I}$. The table measures both properties.

A reconstruction close to $\mathbf{A}$ does not make a nonorthogonal $\mathbf{Q}$ a valid orthogonal factor. Compare the actual orthogonality errors with the reconstruction errors. Householder QR is the library computation; the two Gram–Schmidt routines are printed in full so you can inspect their projection order.

Activity. Choose which of the two reported errors is needed before using $\mathbf{Q}^\mathsf{T}\mathbf{Q}=\mathbf{I}$ in a derivation, and explain your reasoning. Then sweep eps from 10^-1 to 10^-10 and plot both errors for each method. If a null-space method assumes an orthonormal basis, identify the equation that becomes inaccurate when orthogonality is lost.

D4. Symmetric indefinite: the inertia comes back wrong¶

A $4\times4$ KKT matrix, with a positive definite Hessian block $\mathbf{H}=2\mathbf{I}$ and two constraint rows that nearly duplicate each other:

$\mathbf{A}$ is nonsingular for every $\varepsilon\neq0$, so the inertia of $\mathbf{K}$ - the counts of positive, negative and zero eigenvalues - is $(2,2,0)$ exactly, on every row below. The exact column is not an assertion typed in by hand: it is recomputed symbolically with sympy, over the rationals, from the actual double precision entries of $\mathbf{K}$.

The wrong answers printed below are the demonstration. They are what this cell exists to show.

No solver, no barrier parameter, no iteration - just a matrix, a factorization, and a wrong answer about it.

Find the row where the two computed columns disagree with each other, so at least one of them is wrong, and use the exact column to see which. On the machine the handout was written on that is $\varepsilon=10^{-9}$, where ldl is wrong and eigvalsh is right; on yours it may be a different row, or there may be more than one. That variability is itself the lesson: two routines in the same library, asked the same question about the same matrix, return different answers, and which one is wrong where depends on the LAPACK build underneath.

Activity

• An interior point solver reads the inertia of the KKT matrix at every iteration to decide whether to correct the Hessian. What does this table do to that sentence?

• A solver that reads $(3,1,0)$ here concludes the reduced Hessian has the wrong curvature and adds a correction it does not need. Describe the mirror case, where the measurement error points the other way. Which is worse?

• Where did the ill conditioning come from? Note that it is not a scaling accident: it is two constraints that nearly say the same thing, which approaches a constraint-qualification failure as the rows become identical. Each represented matrix in this table still has independent rows: near dependence is not an exact LICQ failure.

Where each of these comes back. None of the four is a curiosity; each is planted here and paid off later, inside an algorithm that has to cope with it.

The algorithmic halves already exist as runnable code: Quasi-Newton Methods for Unconstrained Optimization implements the SR1 update and shows a Cholesky factorization in SciPy, leaving the BFGS update built on it as homework, Inertia-Corrected Newton Method for Equality Constrained NLPs runs inertia correction inside the iteration, and NLP Diagnostics with Degeneracy Hunter is the diagnostic tooling for the last row.

Worked LU, QR, and Cholesky factorizations¶

LU: two triangular solves¶

For $\mathbf{A}=\left[\begin{smallmatrix}4&3\\6&3\end{smallmatrix}\right]$ and $\mathbf{b}=(1,0)^\mathsf{T}$, SciPy’s factors are

First solve $\mathbf{Ly}=\mathbf{P}^\mathsf{T}\mathbf{b}=(0,1)^\mathsf{T}$. Forward substitution gives $y_1=0$ and $\tfrac23y_1+y_2=1$, hence $\mathbf{y}=(0,1)^\mathsf{T}$. Back substitution in $\mathbf{Ux}=\mathbf{y}$ gives $x_2=1$, $6x_1+3x_2=0$, hence $\mathbf{x}=(-1/2,1)^\mathsf{T}$.

QR: projections and a first-row check¶

For $\mathbf{A}\in\mathbb{R}^{m\times n}$, $m\ge n$, with independent columns, the thin QR factorization is $\mathbf{A}=\mathbf{QR}$ with $\mathbf{Q}\in\mathbb{R}^{m\times n}$, $\mathbf{Q}^\mathsf{T}\mathbf{Q}=\mathbf{I}_n$, and nonsingular upper-triangular $\mathbf{R}\in\mathbb{R}^{n\times n}$. For square $\mathbf{A}$, solve $\mathbf{Rx}=\mathbf{Q}^\mathsf{T}\mathbf{b}$. For a tall matrix, the same equation solves the full-column-rank least-squares problem, since the residual decomposes into components in the range of $\mathbf{Q}$ and its orthogonal complement.

Take $\mathbf{A}=\left[\begin{smallmatrix}3&1\\4&2\end{smallmatrix}\right]=[\mathbf{a}_1\ \mathbf{a}_2]$ and $\mathbf{b}=(4,6)^\mathsf{T}$.

1. Normalize the first column: $r_{11}=\|\mathbf{a}_1\|_2=5$, $\mathbf{q}_1=(3/5,4/5)^\mathsf{T}$.

2. Project the second column: $r_{12}=\mathbf{q}_1^\mathsf{T}\mathbf{a}_2=11/5$. Subtract this component: $\mathbf{v}=\mathbf{a}_2-r_{12}\mathbf{q}_1=(-8/25,6/25)^\mathsf{T}$.

3. Normalize the orthogonal remainder: $r_{22}=\|\mathbf{v}\|_2=2/5$, $\mathbf{q}_2=(-4/5,3/5)^\mathsf{T}$. Orthogonality follows from $\mathbf{q}_1^\mathsf{T}\mathbf{v}=\mathbf{q}_1^\mathsf{T}\mathbf{a}_2-r_{12}=0$.

4. Assemble $\mathbf{Q}=\left[\begin{smallmatrix}3/5&-4/5\\4/5&3/5\end{smallmatrix}\right]$ and $\mathbf{R}=\left[\begin{smallmatrix}5&11/5\\0&2/5\end{smallmatrix}\right]$. The dot product of the two columns of $\mathbf{Q}$ is zero.

5. Check the first row of $\mathbf{QR}$: $((3/5)5,\ (3/5)(11/5)-(4/5)(2/5))=(3,1)$, the first row of $\mathbf{A}$.

6. Compute $\mathbf{Q}^\mathsf{T}\mathbf{b}=(36/5,2/5)^\mathsf{T}$. Back substitution gives $x_2=1$, $x_1=1$.

The decomposition $\mathbf{a}_2=r_{12}\mathbf{q}_1+r_{22}\mathbf{q}_2$ is the geometric mechanism: a projection along the first unit direction plus a perpendicular remainder. This two-column Gram–Schmidt calculation explains the factors. Library QR uses Householder transformations; sign choices can change columns of $\mathbf{Q}$ and matching rows of $\mathbf{R}$ while leaving the product unchanged. See Nocedal & Wright (2006), Appendix A, pp.608–609, and Higham (2002), Chapter 19 for numerical analysis.