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.

# Load required Python libraries.
import numpy as np
from scipy import linalg

# Seed the random number generator so this notebook is reproducible
# (Pyomo style guide, section 8).
np.random.seed(0)

Notation used in the references¶

Scalar constants and variables - Greek letters
Vectors - Lower case Roman letters
Matrices - Upper case Roman letters

The notation for lectures will (hopefully) be obvious from context. I sometimes use Greek letters for constants and Roman letters for variables... unless the Greek letters have engineering or scientific meaning.

Matrix Operations¶

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

C_{ij}=\sum_{k=1}^{m}A_{ik}B_{kj}.

(1)

The inner dimensions must match. The transpose interchanges rows and columns, $(A^\mathsf{T})_{ij}=A_{ji}$ . A square matrix is symmetric when $A=A^\mathsf{T}$ , diagonal when $A_{ij}=0$ for $i\ne j$ , and the identity has $I_{ij}=1$ for $i=j$ and zero otherwise.

# Generate random matrices
n = 2
m = 3

Ao = np.random.rand(n, m)
Bo = np.random.rand(m, n)

print("Ao = ")
print(Ao)
print("\nBo = ")
print(Bo)


# What is the dimension of A*B? Try to answer this without the computer.

# Calculate A*B using Python
print("\nAo*Bo =")
print(Ao.dot(Bo))

Ao = 
[[0.5488135  0.71518937 0.60276338]
 [0.54488318 0.4236548  0.64589411]]

Bo = 
[[0.43758721 0.891773  ]
 [0.96366276 0.38344152]
 [0.79172504 0.52889492]]

Ao*Bo =
[[1.40657799 1.08244885]
 [1.15806481 0.98996907]]

Activity

Calculate $A_o^{T}$
Create a square matrix and check if it is symmetric
Create a 3x3 identity matrix

# Transpose
print("transpose(Ao) = \n", Ao.transpose(), "\n")

# Create a square matrix. Is it symmetric?

# Identify matrix
print("I = \n", np.identity(3), "\n")

transpose(Ao) = 
 [[0.5488135  0.54488318]
 [0.71518937 0.4236548 ]
 [0.60276338 0.64589411]] 

I = 
 [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

Determinant¶

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

\det(A)=\sum_{j=1}^{n}(-1)^{i+j}A_{ij}\det(A_{ij}),

(2)

where $A_{ij}$ is the submatrix formed by deleting row $i$ and column $j$ . Useful properties are

\det(AB)=\det(A)\det(B),\qquad \det(A^\mathsf{T})=\det(A),\qquad \det(\alpha A)=\alpha^n\det(A),\qquad \det(I)=1.

(3)

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(A)=0$ .

# Generate random matrices
nd = 2

Ad = np.random.rand(nd, nd)
Bd = np.random.rand(nd, nd)

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

print("Ad = ")
print(Ad)
print("\nBd = ")
print(Bd)
print("\nCd = ")
print(Cd)

Ad = 
[[0.56804456 0.92559664]
 [0.07103606 0.0871293 ]]

Bd = 
[[0.0202184  0.83261985]
 [0.77815675 0.87001215]]

Cd = 
[[ 1  0  0]
 [-2  0  0]
 [ 0  1  1]]

Activity

Verify the properties listed under Some properties for the determinant include
Are $A_d$ , $B_d$ or $C_d$ singular?

# Property 1
print("det(A*B) = ", linalg.det(Ad.dot(Bd)))

print("det(A)*det(B) = ", linalg.det(Ad) * linalg.det(Bd))

det(A*B) =  0.010247347495272154
det(A)*det(B) =  0.010247347495272162

# Property 2

# Property 3

# Property 4

Another Example (for at home)

Expand along the first row of

A=\begin{bmatrix}1&5&0\\2&4&-1\\0&-2&0\end{bmatrix}.

(4)

Then

\det(A)=1\det\begin{bmatrix}4&-1\\-2&0\end{bmatrix}-5\det\begin{bmatrix}2&-1\\0&0\end{bmatrix}+0\det\begin{bmatrix}2&4\\0&-2\end{bmatrix}=-2.

(5)

This is the worked example previously shown as a crop from Lay; it is written natively here so the signs and arithmetic remain searchable and accessible.

Rank¶

The rank of $A\in\mathbb{R}^{n\times m}$ is the dimension of its column space (equivalently, its row space). It is also the size of the largest square submatrix with nonzero determinant. Full rank means $\operatorname{rank}(A)=\min(n,m)$ ; otherwise $A$ is rank deficient.

Activity

Predict the rank of $A_o^T$ , $B_o^T$ , $A_d$ , $B_d$ , $C_d$
Calculate np.linalg.matrix_rank(). Were you correct?

print("rank(A_o) = ", np.linalg.matrix_rank(Ao), "\n")

# Fill in remainder here.

rank(A_o) =  2

Inverse¶

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

Activity

Calculate the inverse of the square matrix $A_d$ and verify that $A_d^{-1} A_d = I$ . (The rectangular matrix $A_o$ has no two-sided inverse.)
Verify that $\det(A_d^{-1}) = 1/\det(A_d)$
Is $Q$ an orthogonal matrix?

# Task 1

# Task 2

# Task 3
# Is this an orthogonal matrix? Yes, no, or need more information?
theta = 0.2
Q = np.array([(np.cos(theta), -np.sin(theta)), (np.sin(theta), np.cos(theta))])

print("Q = ")
print(Q)

Q = 
[[ 0.98006658 -0.19866933]
 [ 0.19866933  0.98006658]]

Solving Linear Systems¶

Consider the linear system $A_l x = b_l$

Al = np.array([(4, 3), (6, 3)])
bl = np.array([1, 0])
print("Al = \n", Al)
print("\nbl = \n", bl)

Al = 
 [[4 3]
 [6 3]]

bl = 
 [1 0]

Explicit Inverse¶

Calculate $x$ by explicitly using $A_l^{-1}$ . Hint: Use linalg.inv().

Ainv = linalg.inv(Al)
print("Ainv = \n", Ainv, "\n")

# Now calculate and print x

Ainv = 
 [[-0.5         0.5       ]
 [ 1.         -0.66666667]]

LU Decomposition¶

Perform LU decomposition on $A_l$ . What structures do the $P$ , $L$ , and $U$ matrices have?

# Create test linear systems
# The transpose is deliberate: it gives a permutation matrix P that is NOT
# symmetric, so the MATLAB and SciPy conventions below give different answers.
A = np.random.rand(4, 4).T
print("A = \n", A, "\n")

b = np.random.rand(4, 1)
print("b = \n", b, "\n")

A = 
 [[0.97861834 0.11827443 0.52184832 0.45615033]
 [0.79915856 0.63992102 0.41466194 0.56843395]
 [0.46147936 0.14335329 0.26455561 0.0187898 ]
 [0.78052918 0.94466892 0.77423369 0.6176355 ]] 

b = 
 [[0.61209572]
 [0.616934  ]
 [0.94374808]
 [0.6818203 ]]

# Perform LU decomposition
P, L, U = linalg.lu(A)

# Permutation matrix
print("P = \n", P)

# Lower triangular matrix
print("L = \n", L)

# Upper triangular matrix
print("U = \n", U)

# Verify result
print("P*L*U = \n", P.dot(L.dot(U)), "\n")

P = 
 [[1. 0. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]]
L = 
 [[1.         0.         0.         0.        ]
 [0.79758282 1.         0.         0.        ]
 [0.81661924 0.63896662 1.         0.        ]
 [0.47156214 0.10299414 0.07659477 1.        ]]
U = 
 [[ 0.97861834  0.11827443  0.52184832  0.45615033]
 [ 0.          0.85033527  0.35801643  0.25381783]
 [ 0.          0.         -0.24024999  0.03375169]
 [ 0.          0.          0.         -0.22504038]]
P*L*U = 
 [[0.97861834 0.11827443 0.52184832 0.45615033]
 [0.79915856 0.63992102 0.41466194 0.56843395]
 [0.46147936 0.14335329 0.26455561 0.0187898 ]
 [0.78052918 0.94466892 0.77423369 0.6176355 ]]

Is P orthogonal?¶

print("inv(P) = \n", linalg.inv(P), "\n")

print("transpose(P) = \n", P.T, "\n")

print("P.T*P = \n", np.matmul(P.T, P), "\n")

inv(P) = 
 [[ 1.  0. -0.  0.]
 [ 0.  0. -0.  1.]
 [ 0.  1. -0.  0.]
 [ 0.  0.  1.  0.]] 

transpose(P) = 
 [[1. 0. 0. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]] 

P.T*P = 
 [[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

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

MATLAB¶

Defines LU decomposition as follows:

P \cdot A = L \cdot U

(6)

Consider

A \cdot x=b

(7)

Using MATLAB LU definition

P\cdot A = L\cdot U

(8)

P^T \cdot P \cdot A = P^T \cdot L \cdot U

(9)

A = P^T \cdot L \cdot U

(10)

Substitute into linear system:

P^T \cdot L \cdot U \cdot x = b

(11)

L \cdot U \cdot x = P \cdot b

(12)

Let $y = U \cdot x$ and substitute.

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

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

Pb = P.dot(b)
print("P*b = ", Pb, "\n")

yLU = linalg.solve(L, Pb)
print("yLU = ", yLU, "\n")

xLU = linalg.solve(U, yLU)
print("xLU = ", xLU, "\n")

# Check: does this x actually solve A*x = b?
print("||A*x - b|| = ", np.linalg.norm(A.dot(xLU) - b), "\n")

P*b =  [[0.61209572]
 [0.94374808]
 [0.6818203 ]
 [0.616934  ]] 

yLU =  [[ 0.61209572]
 [ 0.45555105]
 [-0.10911075]
 [ 0.28973105]] 

xLU =  [[ 0.98256068]
 [ 0.80496658]
 [ 0.27328499]
 [-1.28746252]] 

||A*x - b|| =  0.42382052546128374

SciPy¶

Defines LU decomposition as:

A = P \cdot L \cdot U

(13)

Consider

A \cdot x=b

(14)

SciPy LU definition

A = P \cdot L\cdot U

(15)

Substitute into linear system:

P \cdot L \cdot U \cdot x = b

(16)

P^T \cdot P \cdot L \cdot U \cdot x = P^T \cdot b

(17)

L \cdot U \cdot x = P^T \cdot b

(18)

Let $y = U \cdot x$ and substitute.

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

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

## LU decomposition algorithm.
Pb = (P.T).dot(b)
print("P.T*b = ", Pb, "\n")

yLU = linalg.solve(L, Pb)
print("yLU = ", yLU, "\n")

xLU = linalg.solve(U, yLU)
print("xLU = ", xLU, "\n")

# Check: does this x actually solve A*x = b?
print("||A*x - b|| = ", np.linalg.norm(A.dot(xLU) - b), "\n")

P.T*b =  [[0.61209572]
 [0.6818203 ]
 [0.616934  ]
 [0.94374808]] 

yLU =  [[ 0.61209572]
 [ 0.19362327]
 [-0.00663395]
 [ 0.63567297]] 

xLU =  [[ 1.99078818]
 [ 1.22630486]
 [-0.3692182 ]
 [-2.82470629]] 

||A*x - b|| =  3.6821932062951477e-16

What is the difference? A transpose.

Verify our answer with `linalg.solve`¶

Solve the linear system using linalg.solve

x = linalg.solve(Al, bl)
print("x = \n", x)

x = 
 [-0.5  1. ]

Invertible Matrix Theorem¶

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

$A$ is invertible; its row-reduced form is $I$ ; and it has $n$ pivots.
$Ax=0$ has only the trivial solution, so $\operatorname{Nul}(A)=\{0\}$ and $\dim\operatorname{Nul}(A)=0$ .
The columns of $A$ are linearly independent, form a basis of $\mathbb{R}^n$ , and span $\mathbb{R}^n$ ; equivalently $\operatorname{rank}(A)=n$ .
The map $x\mapsto Ax$ is both one-to-one and onto; therefore $Ax=b$ has a solution for every $b\in\mathbb{R}^n$ .
A left or right inverse exists; $A^\mathsf{T}$ is invertible; $\det(A)\ne0$ ; and zero is not an eigenvalue.
The row space is $\mathbb{R}^n$ , $(\operatorname{Col}A)^\perp=\{0\}$ , $(\operatorname{Nul}A)^\perp=\mathbb{R}^n$ , and $A$ has $n$ nonzero singular values.

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

Eigenvectors and Eigenvalues¶

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

If $A$ is real and symmetric, the spectral theorem gives $A=V\Lambda 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 $x^\mathsf{T}Ax$ . That form is governed by the symmetric part $(A+A^\mathsf{T})/2$ . For rectangular matrices, singular values $\sigma_i\ge0$ —the square roots of the eigenvalues of $A^\mathsf{T}A$ —play the analogous role.

Level sets of a quadratic form, as one picture¶

$Level sets of the quadratic form q(x)=\tfrac12 x^\mathsf{T} A x for symmetric positive definite A, with the eigenvectors drawn as the principal axes. Both panels use the same eigenvectors and the same level q=2; only the eigenvalues change.$

Eigenvalues and eigenvectors of a symmetric matrix have a picture: they are the principal axes of the level sets of $q(x) = \tfrac12 x^\mathsf{T} A x$ . In the rotated coordinates $z = V^\mathsf{T} 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. Note the axis ratio is $\sqrt{\kappa(A)}$ , not $\kappa(A)$ ; by the time you reach the $\kappa \approx 10^{14}$ you will meet in the Condition Number section below, 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.

Ae = np.array([(0, -4, -6), (-1, 0, -3), (1, 2, 5)])
print("Ae = \n", Ae)

Ae = 
 [[ 0 -4 -6]
 [-1  0 -3]
 [ 1  2  5]]

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 $x^\mathsf{T}A_e x>0$ for every nonzero $x$ ?
Form the symmetric part $H=(A_e+A_e^\mathsf{T})/2$ and use its eigenvalues to classify the quadratic form.
Calculate $\det(A_e)$ using only the eigenvalues of $A_e$ .

# Expect negative eigenvalues?
# Yes / No / Cannot tell from inspection?
# Write a sentence to justify your answer.

### Calculate eigenvalues

# Matrix Ae is nonsymmetric
print("Matrix = \n", Ae, "\n")
l, v = linalg.eig(Ae)
print("Eigenvalues = ", l, "\n")
print("Eigenvectors = \n", v, "\n")

# The symmetric part, not Ae's eigenvalues, classifies x.T @ Ae @ x
H = 0.5 * (Ae + Ae.T)
print("Eigenvalues of symmetric part = ", linalg.eigvalsh(H))

Matrix = 
 [[4 3]
 [6 3]] 

Eigenvalues =  [ 7.77200187+0.j -0.77200187+0.j] 

Eigenvectors = 
 [[ 0.62246561 -0.53222953]
 [ 0.78264715  0.8466001 ]]

Singular Value Decomposition¶

Notes will be given in class.

Activity

Calculate singular value decomposition of $A_l$ and $C_d$
What is the rank of each matrix?

# Matrix Al
print("Matrix = \n", Al, "\n")
U, s, Vh = linalg.svd(Al)
print("U = \n", U, "\n")
print("S = \n", s, "\n")
print("Vh = \n", Vh, "\n")

## Rank (from inspecting singular values)?

Matrix = 
 [[4 3]
 [6 3]] 

U = 
 [[-0.59581566 -0.80312122]
 [-0.80312122  0.59581566]] 

S = 
 [8.33557912 0.71980602] 

Vh = 
 [[-0.86400595 -0.50348159]
 [ 0.50348159 -0.86400595]]

Vector and Matrix Norms¶

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

||x|| = 0 \Rightarrow x = 0, \qquad ||\alpha x|| = |\alpha| \cdot ||x||, \qquad ||x + y|| \leq ||x|| + ||y||

(19)

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

$||x||_1 = \sum_{i=1}^n |x_i|$ --- the 1-norm
$||x||_2 = \left( \sum_{i=1}^n x_i^2 \right)^{1/2} = \sqrt{x^T x}$ --- the Euclidean norm
$||x||_{\infty} = \max_i |x_i|$ --- the max norm

Also $|x^T y| \leq ||x||_p \cdot ||y||_q$ whenever $1/p + 1/q = 1$ ; the case $p = q = 2$ is the Cauchy-Schwarz inequality.

Matrix norms are induced by vector norms:

||A|| = \max_{x \neq 0} \frac{||Ax||}{||x||} \qquad \Rightarrow \qquad ||Ax|| \leq ||A|| \cdot ||x||

(20)

The rules for the three induced norms are

$||A||_1 = \max_j \sum_i |A_{ij}|$ --- largest column sum
$||A||_{\infty} = \max_i \sum_j |A_{ij}|$ --- largest row sum
$||A||_2 =$ the largest singular value of $A$

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

||AB|| \leq ||A|| \cdot ||B||

(21)

for any two compatible matrices $A$ and $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).

x = np.random.rand(3, 1)
print("x = \n", x)

y = np.random.rand(3, 1)
print("y = \n", y)

x = 
 [[0.3595079 ]
 [0.43703195]
 [0.6976312 ]]
y = 
 [[0.06022547]
 [0.66676672]
 [0.67063787]]

Activity

Verify the vector norm properties hold for the p-2 norm using $x$ and $y$ . Hint: linalg.norm()
Calculate $||A_l||_1$ , $||A_l||_{\infty}$ and $||A_l||_2$ using both the rules given above and SciPy

Condition Number¶

Activity

Estimate the condition number of $A_l$ by inspecting the SVD results.
Calculate the condition number of $A_l$ using np.linalg.cond()

Notation used in the references¶

Matrix Operations¶

Determinant¶

Rank¶

Inverse¶

Solving Linear Systems¶

Explicit Inverse¶

LU Decomposition¶

Is P orthogonal?¶

MATLAB¶

SciPy¶

Verify our answer with linalg.solve¶

Invertible Matrix Theorem¶

Eigenvectors and Eigenvalues¶

Level sets of a quadratic form, as one picture¶

Singular Value Decomposition¶

Vector and Matrix Norms¶

Condition Number¶

Verify our answer with `linalg.solve`¶