CODEWITHSUNDEEP
pythonnumpylinear-algebra
Carlos Giudice
Carlos Giudice

Reputation: 23

Matrix QR factorization algorithms

I've been trying to visualize QR decomposition in a step by step fashion, but I'm not getting expected results. I'm new to numpy so it'd be nice if any expert eye could spot what I might be missing:

import numpy as np
from scipy import linalg

A = np.array([[12, -51, 4],  
[6, 167, -68],
[-4, 24, -41]])

#Givens

v = np.array([12, 6])

vnorm = np.linalg.norm(v)

W_12 = np.array([[v[0]/vnorm, v[1]/vnorm, 0],
               [-v[1]/vnorm,  v[0]/vnorm, 0],
               [0, 0, 1]])

W_12 * A #this should return a matrix such that [1,0] = 0

#gram-schmidt

A[:,0]

v =  np.linalg.norm(A[:,0]) * np.array([1, 0, 0])

u = (A[:,0] - v) 

u = u / np.linalg.norm(u)

W1 = np.eye(3) - 2 * np.outer(u, u.transpose())

W1 * A #this matrix's first column should look like [a, 0, 0]

any help clarifying the fact that this intermediate results don't show the properties that they are supposed to will be greatly received

Upvotes: 2

Views: 799

Answers (1)

bnaecker
bnaecker

Reputation: 6430

NumPy is designed to work with homogeneous multi-dimensional arrays, it is not specifically a linear algebra package. So by design, the * operator is element-wise multiplication, not the matrix product.

If you want to get the matrix product, there are a few ways:

  1. You can create np.matrix objects, rather than np.ndarray objects, for which the * operator is the matrix product.

  2. You can also use the @ operator, as in W_12 @ A, which is the matrix product.

  3. Or you can use np.dot(W_12, A) or W_12.dot(A), which computes the dot product.

Any one of these, using the data you give, returns the following for Givens rotation:

>>> np.dot(W_12 A)[1, 0]
-2.2204460492503131e-16

And this for the Gram-Schmidt step:

>>> (W1.dot(A))[:, 0]
array([  1.40000000e+01,  -4.44089210e-16,   4.44089210e-16])

Upvotes: 1

Related Questions