Matrix norms, unitary matrices

import numpy as np

Useful tools using numpy

Commented - equal jax operators for advanced users

n, m = 3, 3

# random matrix

np.random.normal(size=(n,m))

# Identity matrix

np.identity(n)

# defining your own matrix

np.array([[4,  8,  15],
           [16, 23, 42],
           [0,  0,  0]])

# matrix multiplication

a = np.random.normal(size=(n,m))
b = np.random.normal(size=(n,m))
x = np.random.normal(size=(n,))
print(np.matmul(a, b))
print()
print(np.matmul(a, x))

a = np.random.normal(size=(n,m))
for i,j in enumerate(a):
    print(i, j[0])

a[:,0]

Matrix Norms Calculation

Important case of operator norms are matrix p-norms, which are defined for \|\cdot\|_* = \|\cdot\|_{**} = \|\cdot\|_p

Among all p-norms three norms are the most common ones:

p = 1, \quad \Vert A \Vert_{1} = \displaystyle{\max_j \sum_{i=1}^n} |a_{ij}|

p = 2, \quad spectral norm, denoted by \Vert A \Vert_2

p = \infty, \quad \Vert A \Vert_{\infty} = \displaystyle{\max_i \sum_{j=1}^m} |a_{ij}|

Let us work with = 1 and inf now.

enumerate(a)

Lets write our own functions to calculate 0 and inf norms:

def p_1_norm(a):
    n,m = a.shape
    max = 0
    for i in range(m):
        max = np.max([np.sum(np.abs(a[:,[i]])), max])
    return max

def p_inf_norm(a):
    n,m = a.shape
    max = 0
    for i in range(n):
        max = np.max([np.sum(np.abs(a[i])), max])
    return max

m, n = 3, 3
a = np.random.normal(size=(n,m)) #Random n x m matrix
print(a)
print("p_1 norm: ", p_1_norm(a))
print("p_inf norm: ", p_inf_norm(a))

Now, make up a toy example matrix for which these two are not equal:

b = np.array([[1, 0, 0],
               [0, 1, 0],
               [0, 0, 1]])

print("p_1 norm: ", p_1_norm(b))
print("p_inf norm: ", p_inf_norm(b))
print(p_1_norm(b) == p_inf_norm(b))

Creative task: invent your own matrix norm and check that it satisfies three norm requirements

\Vert \cdot \Vert is called a matrix norm if it is a vector norm on the vector space of n \times m matrices and:

\|A\| \geq 0 and if \|A\| = 0 then A = O

\|\alpha A\| = |\alpha| \|A\|

\|A+B\| \leq \|A\| + \|B\| (triangle inequality)

Additionally some norms can satisfy the submultiplicative property: \Vert A B \Vert \leq \Vert A \Vert \Vert B \Vert

def super_norm(a):
    
    n,m = a.shape
    sum = 0
    for i in range(n):
        sum += np.abs(a[i,i])
    
    return sum

n, m = 50, 50

Now, lets check all requirements. Attention! If some conditions seem to not work, recall previous lecture to understand and fix the problem

is_nonnegative = True

for i in range(42):
    a = np.random.normal(size=(n,m))
    
    if super_norm(a) < 0:
        is_nonnegative = False
        
print("is_nonnegative: ", is_nonnegative)

respects_zero = True

a = np.zeros((n,m))

if super_norm(a) != 0:
    respects_zero = False
    
print("respects_zero: ", respects_zero)

respects_scalar_mult = True

for i in range(42):
    a = np.random.normal(size=(n,m))
    alpha = np.abs(np.random.normal())
    
    if np.abs(super_norm(alpha * a) - alpha * super_norm(a)) > 0.001:
        respects_scalar_mult = False
        
print("respects_scalar_mult: ", respects_scalar_mult)

triangle_inequality = True

for i in range(42):
    a = np.random.normal(size=(n,m))
    b = np.random.normal(size=(n,m))
    
    if super_norm(a + b) > super_norm(a) + super_norm(b):
        triangle_inequality = False
        
print("triangle_inequality: ", triangle_inequality)

# Not required to be a norm

is_submultiplicative = True

for i in range(42):
    a = np.random.normal(size=(n,m))
    b = np.random.normal(size=(n,m))
    
    if super_norm(np.matmul(a, b)) > super_norm(a) * super_norm(b):
        is_submultiplicative = False
        
print("is_submultiplicative: ", is_submultiplicative)

Householder Martices

Householder matrix is the matrix of the form: H \equiv H(v) = I - 2 vv^*, where v is an n \times 1 column and v^* v = 1.

It is also a reflection: Hx = x - 2(v^* x) v

Attention! If it does not work, remember about vector norm

Build your own from a vector

def build_householder(v):
    # v - vector of size n
    v = v/np.linalg.norm(v)
    a = np.identity(v.size)
    
    H = a - 2*np.outer(v, v)
    
    return H

v = np.random.normal(size=(3))
print(v)
h = build_householder(v)
print(h)

See how it reflects vectors

v = np.random.normal(size=(3))
h = build_householder(v)

x = np.random.normal(size=(3))
print(x)
print(np.matmul(h, x))

v = np.array([0,  1,  -1])
h = build_householder(v)
x = np.array([0,  1,  0])
print(np.round(h , decimals=2))
print(x)
print(np.round(h @ x.T, decimals=2))

Optional: check that it indeed is also a reflection

Hx = x - 2(v^* x) v

v = np.random.normal(size=(3))
v = v/np.linalg.norm(v)
h = build_householder(v)

x = np.random.normal(size=(3))

v = np.array([0,  1,  -1])
v = v/np.linalg.norm(v)
h = build_householder(v)
x = np.array([0,  1,  0])

hx = h @ x.T

reflected =  x - 2 * np.dot(v.T,x) * v

print("initial vector: ", x)
print("transofrmed by matrix: ", hx)
print("reflected by vector: ", reflected)

Check unitarity

Use numpy tools to check if the matrix is unitary using formula U* x U = I

n = 3

v = np.random.normal(size=(n))
h = build_householder(v)

print(np.round(h @ h, decimals=2))

def householder_transform(A):
    """
    Transforms the matrix A into an upper triangular matrix using Householder reflections.
    
    Parameters:
        A (numpy.ndarray): The matrix to be transformed.
    
    Returns:
        R (numpy.ndarray): The upper triangular matrix after applying Householder transformations.
    """
    A = A.copy()
    m, n = A.shape
    for j in range(min(m, n)):
        # Create the vector x for the current column
        x = A[j:, j]
        
        # Calculate the norm of x and the Householder vector v
        norm_x = np.linalg.norm(x)
        if norm_x == 0:
            continue
        sign = -1 if x[0] < 0 else 1
        v = x.copy()
        v[0] += sign * norm_x  # Adjust the first element of v for the reflection
        v /= np.linalg.norm(v)  # Normalize v
        
        # Apply the Householder transformation to A[j:, j:]
        A[j:, j:] -= 2 * np.outer(v, v @ A[j:, j:])
    
    return A

# Example matrix
A = np.array([
    [4, 1, -2, 2],
    [1, 2, 0, 1],
    [-2, 0, 3, -2],
    [2, 1, -2, -1]
], dtype=float)

R = householder_transform(A)
print("Upper triangular matrix R:\n", R)

Bonus task: check that it also preserves the norm. You can check it for your own custom norm if you created one!

\frac{\Vert x - \widehat{x} \Vert}{\Vert x \Vert} \leq \varepsilon.

\frac{\Vert y - \widehat{y} \Vert}{\Vert y \Vert } = \frac{\Vert U ( x - \widehat{x}) \Vert}{\Vert U x\Vert} \leq \varepsilon.

v = np.random.normal(size=(n))
v = v/np.linalg.norm(v)
h = build_householder(v)
epsilon = 0.001

x = np.random.normal(size=(n))
x_hat = x + epsilon * np.random.normal(size=(n)) # approximaton of x
y = x - x_hat                                    # error of approximation

initial_error = np.linalg.norm(y)/np.linalg.norm(x)
transformed_error = np.linalg.norm(h @ y.T)/np.linalg.norm(h @ x.T)
        
print("initial error:     ", initial_error)
print("transformed error: ", transformed_error)