STA 243 Computational Statistics

Discussion 2: Linear Algebra in Python

TA: Tesi Xiao

In Python, we usually use NumPy to implement the matrix computations for the sake of efficiency. The package NumPy provides several powerful classes and methods for numerical computations.

To be specific, the class np.ndarray (np.array) is commonly used for matrix computations. Please check out this manual.

Vector - 1d array

import numpy as np
vec = np.array([1,2,3])
vec.shape
(3,)

(3,) is a tuple with only one element. The 1-d array has only the length attribute.

len(vec)
3

Matrix - 2d array

mat = np.array([[1,2,3],[4,5,9],[7,8,9],[10,11,12]])
mat
array([[ 1,  2,  3],
       [ 4,  5,  9],
       [ 7,  8,  9],
       [10, 11, 12]])
mat.shape
(4, 3)

(4, 3) is a tuple with two elements, which corresponds to # of rows, # of columns.

Tensor - multidimensional array

Tensors are a generalized data structure, which have a wide range of applications in modern machine learning. We will later discuss it in the section of deep neural networks. For example, in PyTorch, we will use tensors to encode the inputs and outputs of a model, as well as the model’s parameters. Then, tensors can run on GPUs or other specialized hardware to accelerate computing.

Matrix Computations

1. Matrix Addition, Substraction, Element-wise Multiplication & Division

A = np.array([[1,2],[3,4]])
B = np.array([[2,3],[4,5]])
A+B
array([[3, 5],
       [7, 9]])
A-B
array([[-1, -1],
       [-1, -1]])
A*B
array([[ 2,  6],
       [12, 20]])
A/B
array([[0.5       , 0.66666667],
       [0.75      , 0.8       ]])

+, -, *, / will do element-wise addition and substraction. If the shapes of two arrays are different, this operation cannot be done with an error message. However, the matrix-scalar calculations are acceptable.

A + 1
array([[2, 3],
       [4, 5]])
A * 2
array([[2, 4],
       [6, 8]])

2. Matrix Muplication

When computing the matrix muplication using NumPy and the operator @, make sure that the shapes of two operands follow the rules, i.e., (m,p) @ (p,n) $\rightarrow$ (m,n) .

  • Matrix-Vector Product (2darray-1darray Product)
b = np.array([1,1,1])
mat @ b
array([ 6, 18, 24, 33])
c = np.array([1,1,1,1])
c @ mat
array([22, 26, 33])

As long as two shapes follow the matrix multiplication rules, it returns a 1d-array.

  • Matrix-Matrix Product (2darray-2darray Product)
b = b.reshape(-1,1)
b, b.shape
(array([[1],
        [1],
        [1]]), (3, 1))
mat @ b
array([[ 6],
       [18],
       [24],
       [33]])
A @ B
array([[10, 13],
       [22, 29]])

It returns a 2d-array following the multiplication rules.

  • Matrix Transpose
mat.T
array([[ 1,  4,  7, 10],
       [ 2,  5,  8, 11],
       [ 3,  9,  9, 12]])

3. Eigenvalues and Eigenvectors

np.linalg.eig

# set a random seed
np.random.seed(2021)
# Create a 3 by 3 random matrix
A = np.random.rand(3,3)

# Compute eigenvalues and eigenvectors
eigenvals, eigenvecs = np.linalg.eig(A)

eigenvals
array([1.46993771, 0.28721889, 0.50822548])
eigenvecs
array([[-0.55986774, -0.76336389, -0.14005302],
       [-0.54059029,  0.42411145, -0.16626493],
       [-0.62794128, -0.48724229,  0.97608459]])

Note that only if all the eigenvectors are linearly independent, we have the eigendecomposition $A = U\Sigma U^\top$.

4. Sigular Value Decomposition (SVD)

np.linalg.svd

$A_{m\times n} = U_{m\times m} S_{m\times n} V^\top_{n\times n}$

# Create a 5 by 3 random matrix
A = np.random.randn(5, 3)
# Singular Value Decomposition
U, S, Vt = np.linalg.svd(A, full_matrices=True)
U, U.shape
(array([[ 0.22073128, -0.03043587,  0.16712615, -0.89000162,  0.36099491],
        [-0.11471826,  0.92761524, -0.33552023, -0.09536938,  0.06856048],
        [ 0.22007922,  0.12998267,  0.32382857,  0.43547724,  0.80010267],
        [ 0.68167637, -0.16826726, -0.70307225,  0.0763847 ,  0.08281451],
        [-0.65192015, -0.30560473, -0.51021515, -0.057678  ,  0.46686145]]),
 (5, 5))
S, S.shape
(array([3.36166521, 2.40725415, 1.48282724]), (3,))

S here is a 1d-array containing non-zero diagonal entries.

Vt, Vt.shape
(array([[-0.5160804 ,  0.28883446, -0.80637193],
        [ 0.19164251,  0.95649989,  0.21995704],
        [ 0.83482583, -0.04101963, -0.5489838 ]]), (3, 3))

5. Norms

np.linalg.norm(x, ord=None, axis=None, keepdims=False)

Below are several useful norms in this class

ordmatrix normvector norm
NoneFrobenius norm2-norm
1max(sum(abs(x), axis=0))1-norm
22-norm (largest singular value)2-norm
infmax(sum(abs(x), axis=1))max(abs(x))
0sum(x != 0)
np.linalg.norm(A) # Frobenius norm
4.392543930507632
np.linalg.norm(b, 1) # L1-norm for vector b 
3.0