Householder qr factorization python

Saved searches

Use saved searches to filter your results more quickly

You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session. You switched accounts on another tab or window. Reload to refresh your session.

Python implementation of QR decomposition using Householder transformations

TayssirDo/QR-decomposition

This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

Sign In Required

Please sign in to use Codespaces.

Launching GitHub Desktop

If nothing happens, download GitHub Desktop and try again.

Launching GitHub Desktop

If nothing happens, download GitHub Desktop and try again.

Launching Xcode

If nothing happens, download Xcode and try again.

Launching Visual Studio Code

Your codespace will open once ready.

There was a problem preparing your codespace, please try again.

Latest commit

Git stats

Files

Failed to load latest commit information.

README.md

We can compute the QR decomposition by Householder transformations, Givens transformations or by Gram Schmidt orthogonalization. Here, we compute Householder QR.

If has full column rank, i.e, Rank( )= then has the QR decomposition:

where is orthogonal and is nonsingular upper triangular.

np.random.seed(20) A = np.random.randn(5, 3) print('----- Matrix A: -----\n'+ str(A)+'\n') #compute the QR decomposition Q, R = qr_householder(A) print('\033[1;31;31m Compute A=QR: \033[1;m \n'+ str(Q @ R)) print(np.allclose(A, Q @ R))

1. G.H. Golub and C.F. Van Loan. Matrix Computations. The Johns Hopkins University Press, Baltimore, Maryland, 4th edition, 2013.

About

Python implementation of QR decomposition using Householder transformations

Источник

Hsankesara / qr_householder.py

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters

«»»

This is the code for QR factorization using Householder Transformation.

This program is made in python 3.5.3 but will be compatible to any python 3.4+ version

We used numpy library for matrix manipulation.

Install numpy using ** pip3 install numpy ** command on terminal.

To run the code write ** python3 qr_householder.py ** on terminal

User has to give dimension of the matrix as input in space separated format and matrix will be generated randomly.

QR factorization can be done for both square and non-square matrices and hence the code supports both.

«»»

import numpy as np

def column_convertor ( x ):

«»»

Converts 1d array to column vector

«»»

x . shape = ( 1 , x . shape [ 0 ])

return x

def get_norm ( x ):

«»»

Returns Norm of vector x

«»»

return np . sqrt ( np . sum ( np . square ( x )))

def householder_transformation ( v ):

«»»

Returns Householder matrix for vector v

«»»

size_of_v = v . shape [ 1 ]

e1 = np . zeros_like ( v )

e1 [ 0 , 0 ] = 1

vector = get_norm ( v ) * e1

if v [ 0 , 0 ] < 0 :

vector = — vector

u = ( v + vector ). astype ( np . float32 )

H = np . identity ( size_of_v ) — (( 2 * np . matmul ( np . transpose ( u ), u )) / np . matmul ( u , np . transpose ( u )))

return H

def qr_step_factorization ( q , r , iter , n ):

«»»

Return Q and R matrices for iter number of iterations.

«»»

v = column_convertor ( r [ iter :, iter ])

Hbar = householder_transformation ( v )

H = np . identity ( n )

H [ iter :, iter :] = Hbar

r = np . matmul ( H , r )

q = np . matmul ( q , H )

return q , r

def main ():

n , m = list ( map ( int , input ( «Write size of the matrix in space separated format \n » ). split ()))

A = np . random . rand ( n , m )

print ( ‘The random matrix is \n ‘ , A )

Q = np . identity ( n )

R = A . astype ( np . float32 )

for i in range ( min ( n , m )):

# For each iteration, H matrix is calculated for (i+1)th row

Q , R = qr_step_factorization ( Q , R , i , n )

min_dim = min ( m , n )

R = np . around ( R , decimals = 6 )

R = R [: min_dim , : min_dim ]

Q = np . around ( Q , decimals = 6 )

print ( ‘A after QR factorization’ )

print ( ‘R matrix’ )

print ( R , ‘ \n ‘ )

print ( ‘Q matrix’ )

print ( Q )

if __name__ == «__main__» :

main ()

Источник

Householder qr factorization python

QR Decomposition with Python and NumPy

QR Decomposition with Python and NumPy

This article will discuss QR Decomposition in Python. In previous articles we have looked at LU Decomposition in Python and Cholesky Decomposition in Python as two alternative matrix decomposition methods. QR Decomposition is widely used in quantitative finance as the basis for the solution of the linear least squares problem, which itself is used for statistical regression analysis.

One of the key benefits of using QR Decomposition over other methods for solving linear least squares is that it is more numerically stable, albeit at the expense of being slower to execute. Hence if you are performing a large quantity of regressions as part of a trading backtest, for instance, you will need to consider very extensively whether QR Decomposition is the best fit (excuse the pun).

For a square matrix $A$ the QR Decomposition converts $A$ into the product of an orthogonal matrix $Q$ (i.e. $Q^TQ=I$) and an upper triangular matrix $R$. Hence:

There are a few different algorithms for calculating the matrices $Q$ and $R$. We will outline the method of Householder Reflections, which is known to be more numerically stable the the alternative Gramm-Schmidt method. I’ve outlined the Householder Reflections method below.

Note, the following explanation is an expansion of the extremely detailed article on QR Decomposition using Householder Reflections over at Wikipedia.

A Householder Reflection is a linear transformation that enables a vector to be reflected through a plane or hyperplane. Essentially, we use this method because we want to create an upper triangular matrix, $R$. The householder reflection is able to carry out this vector reflection such that all but one of the coordinates disappears. The matrix $Q$ will be built up as a sequence of matrix multiplications that eliminate each coordinate in turn, up to the rank of the matrix $A$.

The first step is to create the vector $\mathbb$, which is the $k$-th column of the matrix $A$, for step $k$. We define $\alpha = -sgn(\mathbb_k)(||\mathbb||)$. The norm $||\cdot||$ used here is the Euclidean norm. Given the first column vector of the identity matrix, $I$ of equal size to $A$, $\mathbb_1 = (1,0. 0)^T$, we create the vector $\mathbb$:

\begin \mathbb = \mathbb + \alpha \mathbb_1 \end

Once we have the vector $\mathbb$, we need to convert it to a unit vector, which we denote as $\mathbb$:

Now we form the matrix $Q$ out of the identity matrix $I$ and the vector multiplication of $\mathbb$:

\begin Q = I — 2 \mathbb \mathbb^T \end

$Q$ is now an $m\times m$ Householder matrix, with $Q\mathbb = \left( \alpha, 0, . 0\right)^T$. We will use $Q$ to transform $A$ to upper triangular form, giving us the matrix $R$. We denote $Q$ as $Q_k$ and, since $k=1$ in this first step, we have $Q_1$ as our first Householder matrix. We multiply this with $A$ to give us:

\begin Q_1A = \begin \alpha_1&\star&\dots&\star\\ 0 & & & \\ \vdots & & A’ & \\ 0 & & & \end \end

The whole process is now repeated for the minor matrix $A’$, which will give a second Householder matrix $Q’_2$. Now we have to «pad out» this minor matrix with elements from the identity matrix such that we can consistently multiply the Householder matrices together. Hence, we define $Q_k$ as the block matrix:

Once we have carried out $t$ iterations of this process we have $R$ as an upper triangular matrix:

\begin R = Q_t . Q_2 Q_1 A \end

$Q$ is then fully defined as the multiplication of the transposes of each $Q_k$:

\begin Q = Q^T_1 Q^T_2 . Q^T_t \end

This gives $A=QR$, the QR Decomposition of $A$.

To calculate the QR Decomposition of a matrix $A$ with NumPy/SciPy, we can make use of the built-in linalg library via the linalg.qr function. This is significantly more efficient than using a pure Python implementation:

import pprint import scipy import scipy.linalg # SciPy Linear Algebra Library A = scipy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]) # From the Wikipedia Article on QR Decomposition Q, R = scipy.linalg.qr(A) print "A:" pprint.pprint(A) print "Q:" pprint.pprint(Q) print "R:" pprint.pprint(R) 

The output of the QR decomposition includes $A$, $Q$ and $R$. As a basic sanity check we can see that $R$ is in fact an upper triangular matrix:

A: array([[ 12, -51, 4], [ 6, 167, -68], [ -4, 24, -41]]) Q: array([[ 0.85714286, -0.39428571, -0.33142857], [ 0.42857143, 0.90285714, 0.03428571], [-0.28571429, 0.17142857, -0.94285714]]) R: array([[ 14., 21., -14.], [ -0., 175., -70.], [ 0., -0., 35.]]) 

You aren’t likely to ever need a pure Python implementation of QR Decomposition (homework notwithstanding), but I feel that it is helpful to gain an understanding of the Householder Reflections algorithm, so I have written my own implementation:

from math import sqrt from pprint import pprint def mult_matrix(M, N): """Multiply square matrices of same dimension M and N""" # Converts N into a list of tuples of columns tuple_N = zip(*N) # Nested list comprehension to calculate matrix multiplication return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M] def trans_matrix(M): """Take the transpose of a matrix.""" n = len(M) return [[ M[i][j] for i in range(n)] for j in range(n)] def norm(x): """Return the Euclidean norm of the vector x.""" return sqrt(sum([x_i**2 for x_i in x])) def Q_i(Q_min, i, j, k): """Construct the Q_t matrix by left-top padding the matrix Q with elements from the identity matrix.""" if i < k or j < k: return float(i == j) else: return Q_min[i-k][j-k] def householder(A): """Performs a Householder Reflections based QR Decomposition of the matrix A. The function returns Q, an orthogonal matrix and R, an upper triangular matrix such that A = QR.""" n = len(A) # Set R equal to A, and create Q as a zero matrix of the same size R = A Q = [[0.0] * n for i in xrange(n)] # The Householder procedure for k in range(n-1): # We don't perform the procedure on a 1x1 matrix, so we reduce the index by 1 # Create identity matrix of same size as A I = [[float(i == j) for i in xrange(n)] for j in xrange(n)] # Create the vectors x, e and the scalar alpha # Python does not have a sgn function, so we use cmp instead x = [row[k] for row in R[k:]] e = [row[k] for row in I[k:]] alpha = -cmp(x[0],0) * norm(x) # Using anonymous functions, we create u and v u = map(lambda p,q: p + alpha * q, x, e) norm_u = norm(u) v = map(lambda p: p/norm_u, u) # Create the Q minor matrix Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in xrange(n-k)] for j in xrange(n-k) ] # "Pad out" the Q minor matrix with elements from the identity Q_t = [[ Q_i(Q_min,i,j,k) for i in xrange(n)] for j in xrange(n)] # If this is the first run through, right multiply by A, # else right multiply by Q if k == 0: Q = Q_t R = mult_matrix(Q_t,A) else: Q = mult_matrix(Q_t,Q) R = mult_matrix(Q_t,R) # Since Q is defined as the product of transposes of Q_t, # we need to take the transpose upon returning it return trans_matrix(Q), R A = [[12, -51, 4], [6, 167, -68], [-4, 24, -41]] Q, R = householder(A) print "A:" pprint(A) print "Q:" pprint(Q) print "R:" pprint(R)

The output from the Householder method implemented in pure Python is given below:

A: [[12, -51, 4], [6, 167, -68], [-4, 24, -41]] Q: [[0.8571428571428571, 0.39428571428571435, -0.33142857142857135], [0.4285714285714286, -0.9028571428571429, 0.034285714285714114], [-0.28571428571428575, -0.17142857142857126, -0.942857142857143]] R: [[13.999999999999998, 21.00000000000001, -14.000000000000004], [-5.506706202140776e-16, -175.00000000000003, 70.0], [3.0198066269804245e-16, -3.552713678800501e-14, 35.000000000000014]] 

You can see that we get the same answers as above in the SciPy implementation, albeit with a few more significant figures! One has to be extremely careful in numerical algorithms when dealing with floating point arithmetic, but that is a discussion for another day.

QSAlpha

Join the QSAlpha research platform that helps fill your strategy research pipeline, diversifies your portfolio and improves your risk-adjusted returns for increased profitability.

The Quantcademy

Join the Quantcademy membership portal that caters to the rapidly-growing retail quant trader community and learn how to increase your strategy profitability.

Successful Algorithmic Trading

How to find new trading strategy ideas and objectively assess them for your portfolio using a Python-based backtesting engine.

Advanced Algorithmic Trading

How to implement advanced trading strategies using time series analysis, machine learning and Bayesian statistics with R and Python.

Источник