$A = QR$

Code implementation

We already have the method to give us a projection matrix. So here we take $I - P$ as our projection matrix. After transposing $A$ (to make accessing the columns easier) we loop over each column and make the others columns orthogonal. Then we move onto the next column and repeat:

def qr(self):
    if self.is_singular():
        print('Matrix is singular!')
        return self, None, None

    A = self.tr()
    Q = dc(A)
    I = eye(A.size())
    # projection orthogonal to column
    for col in range(Q.size(0)-1):
        Col = dc(Mat([Q.data[col]]))
        P, _ = Col.tr().projection()
        P = I.subtract(P)
        # project and put into matrix Q
        for col2 in range(col+1, Q.size(0)):
            Col = dc(Mat([Q.data[col2]]))
            q = P.multiply(Col.tr()).tr()
            Q.data[col2] = q.data[0]

The next step is to make the length of each column 1. Which we do by dividing by it's current length:

# normalise to unit length
for x, q in enumerate(Q.data):
    q = Mat([q])
    q = q.norm()
    Q.data[x] = q.data[0]

The last step is to transpose $A$ back to how it was originally, and to generate $R$ by multiplying $Q^TA$. Note that because we worked from $A^T$ from the beginning, our $Q$ is already transposed, so we just multiply this with $A$ and then transpose $Q$ before returning:

    A = A.tr()
    R = Q.multiply(A)
    Q = Q.tr()
    A = Q.multiply(R)

    return A, Q, R

Demo

We create a matrix, call the qr method and print the results. Note that we print the $A$ computed from $QR$ at two levels of rounding precision. We begin to see rounding errors beyond 13 decimal places:

import linalg as la

A = la.Mat([[1, 2, 3],
            [2, 4, 12],
            [2, 8, 0]])

A, Q, R = A.qr()
la.print_mat(Q, 2)
la.print_mat(R, 2)
la.print_mat(A, 13)
la.print_mat(A, 20)

Outputs:

>>> la.print_mat(Q, 2)
[0.33, -0.3, -0.89]
[0.67, -0.6, 0.45]
[0.67, 0.75, 0.0]

>>> la.print_mat(R, 2)
[3.0, 8.67, 9.0]
[0.0, 2.98, -8.05]
[0.0, 0.0, 2.68]

>>> la.print_mat(A, 13)
[1.0, 2.0, 3.0]
[2.0, 4.0, 12.0]
[2.0, 8.0, 0.0]

>>> la.print_mat(A, 20)
[0.9999999999999989, 1.9999999999999973, 2.9999999999999956]
[2.0000000000000004, 4.0, 12.000000000000004]
[2.0, 8.0, 0.0]

< Projection and regression

Eigenvalues and eigenvectors >

back to project main page
back to home

---