Basic linear algebra using numpy

Dot product

Input

import numpy as np
x = np.array([3, 4, 1])
y = np.array([2, -3, 4])
np.dot(x, y)

Output

-2

Angle between two vectors:

$$\cos(\theta) = \frac{\bf{x} \cdot \bf{y}}{|\bf{x}| |\bf{y}|}$$

Input

np.dot(x, y)/(np.linalg.norm(x) * np.linalg.norm(y))

Output

-0.07283570407292297

Matrix multiplication

The inner dimensions of matrices must match for multiplication. The dimension or output matrix is the outer dimensions, e.g., a $2 \times 3$ multiplied with a $3 \times 3$ matrix will result in a $2 \times 3$ matrix.

Matrix multiplication is not commutative (in general): $AB \neq BA$

Input

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

Output

array([[ 1, -2,  1],
       [ 2, -4,  5]])

Input

B = np.array([[3, 4, 2], [1, -2, 0], [2, 1, -3]]); B

Output

array([[ 3,  4,  2],
       [ 1, -2,  0],
       [ 2,  1, -3]])

Input

np.matmul(A, B)

Output

array([[  3,   9,  -1],
       [ 12,  21, -11]])

Transpose

$$(AB)^T = B^T A^T$$

Input

AT = A.transpose(); AT

Output

array([[ 1,  2],
       [-2, -4],
       [ 1,  5]])

Input

(np.matmul(A, B)).transpose()

Output

array([[  3,  12],
       [  9,  21],
       [ -1, -11]])

Input

np.matmul(B.transpose(), A.transpose())

Output

array([[  3,  12],
       [  9,  21],
       [ -1, -11]])

If $u$, $v~\epsilon~\rm{R}^d$; $\bf{u} \cdot \bf{v} = u^Tv$

Input

u = np.array([3, 4, 1]); u

Output

array([3, 4, 1])

Input

v = np.array([1, 5, 3]); v

Output

array([1, 5, 3])

Input

np.matmul(u.transpose(), v)

Output

26

Input

np.dot(u, v)

Output

26

$M$ is symmetric if $M = M^T$ Diagonal matrices are symmetric.

Determinant

For a $2 \times 2$ matrix $A = \begin{pmatrix}a & b\\ c & d \end{pmatrix}$, $det (A) = ad -bc$

For a diagonal matrix, $det(A)$ is the product of diagonal elements.

$$det(AB) = det(A) det(B)$$

Input

A = np.array([[1, 2, 4], [-2, 0, 3], [5, -1, -2]]); A

Output

array([[ 1,  2,  4],
       [-2,  0,  3],
       [ 5, -1, -2]])

Input

np.linalg.det(A)

Output

33.000000000000014

Inverse

$B$ is inverse of $A$ ($A$ must be a square matrix) if $AB = BA = I$

$$B = A^{-1}$$

A matrix for which inverse does not exist is called singular matrix. $A$ is invertible if and only if $|A| \neq 0$.

Input

np.linalg.inv(A)

Output

array([[ 9.09090909e-02, -2.22044605e-17,  1.81818182e-01],
       [ 3.33333333e-01, -6.66666667e-01, -3.33333333e-01],
       [ 6.06060606e-02,  3.33333333e-01,  1.21212121e-01]])

System of linear equations

The central problem of linear algebra is solving the system of linear equations. There are two main methods to solve linear equations: (1) method of elimination and (2) Cramer's rule: method of determinants.

Let us consider the case of $n$ linear equations with $n$ unknowns. In case of elimination, the multiples of first equation is subtracted from the remaining equations to eliminate the first unknown. This leaves a smaller system of $(n-1)$ equations with $(n-1)$ unknowns. The method is repeated until we are left with one equation with one unknown. Then we substitute the solution and find the other unknowns in reversed ordered. This method is generally used in practice to solve system of linear equation, and known as Gaussian elimination.

$$2 x_1 + x_2 + x_3 = 1$$ $$6 x_1 + 2 x_2 + x_3 = -1$$ $$-2 x_1 + 2 x_2 + x_3 = 7$$

Matrix equation: $A\textbf{x} = \textbf{b}$

Solution: $\textbf{x} = A^{-1} \textbf{b}$

$$\begin{bmatrix}2 & 1 & 1 \\ 6 & 2 & 1 \\ -2 & 2 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 \\ -1 \\ 7 \end{bmatrix}$$

Input

A = np.array([[2, 1, 1], [6, 2, 1], [-2, 2, 1]])
b = np.array([1, -1, 7])
x = np.linalg.solve(A, b)
print(x)

Output

[-1.  2.  1.]

LU factorization

After the Gauss elimination, we are left with an upper triangular matrix ($U$).

$$U\textbf{x} = \textbf{c}$$

The matrix operation that relates $A$ and $U$ ($\textbf{b}$ and $\textbf{c}$) is found to be a lower triangular matrix $L$.

$$A = LU$$ $$A\textbf{x} = LU\textbf{x} = L\textbf{c} = \textbf{b}$$

Input

import scipy.linalg as la
p, l, u = la.lu(A)
print("L=", l)
print("U=", u)

Output

L= [[ 1.          0.          0.        ]
 [-0.33333333  1.          0.        ]
 [ 0.33333333  0.125       1.        ]]
U= [[6.         2.         1.        ]
 [0.         2.66666667 1.33333333]
 [0.         0.         0.5       ]]

Input

np.matmul(l, u)

Output

array([[ 6.,  2.,  1.],
       [-2.,  2.,  1.],
       [ 2.,  1.,  1.]])

Notice that the order of the rows changed and this information is contained in the permutaion matrix $P$.

The diagonal elements of $L$ are always 1. There is another form of decomposition $A = LDU$, where both $L$ and $U$ has diagonal elements as 1, and $D$ is the diagonal matrix of pivots. For a nonsingular matrix $A$, the $LU$ and $LDU$ factorization are unique.

Rank of a matrix is the number of pivots or the number of non-zero rows in $U$. It represents the number of independent rows in the matrix. Note that the same number of columns in $A$ are linearly independent as well.

Projections

The projection of a vector $\textbf{b}$ onto $\textbf{a}$ is:

$$p = \frac{a^T b}{a^T a} a$$

The projection matrix $P = \frac{a a^T}{a^T a}$ so that $p = Pb$.

$P$ is a symmetric matrix, and $P^2 = P$

Eigenvalues and eigenvectors

$$A \textbf{x} = \lambda \textbf{x}$$

where the numbers $\lambda$ are the eigenvalue of $A$. $\textbf{x}$ are the special vectors that does not change direction after multiplied by $A$.

• $\det(A - \lambda I) = 0$

• If $A \textbf{x} = \lambda \textbf{x}$ then: $A^2 \textbf{x} = \lambda^2 \textbf{x}$ and $A^{-1} \textbf{x} = \lambda^{-1} \textbf{x}$

• $\det(A) = (\lambda_1) (\lambda_2) \cdots (\lambda_n)$

Resources

• Linear Algebra and Its Applications by Gilbert Strang.
• Matrix Analysis and Applied Linear Algebra by Carl D. Meyer.