## Gramm Schmidt orthogonalization algorithm in Python

I am currently teaching Econometrics and I am encouraging my students to apply open source software like R and Python and Gretl. The sad thing is that they have access to the computer lab only in the third year of their course. Anyway, in the course of writing R and Python programs to illustrate concepts, I have been more appreciative that linear algebra techniques are heavily used behind the least squares method for determining coefficients of a model equation like the QR technique which we will expound further in this blog.

One direction to take in finding the QR decomposition of a matrix A is using the gramm-schmidt orthogonalization procedure to create an orthonormal matrix B from the matrix to A and having the property the dot product of any two column vectors of B is either one or zero and that t(B) B = I_n where the last matrix is the identity matrix.

The algorithm is described in Wikipedia Gram-Schmidt article. Here is a Python implementation which took me few hours of rest periods between teaching.

^{?}View Code PYTHON

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 | """ File gramm.py Author Ernesto P. Adorio, Ph.D. U.P. Clarkfield, Pampanga Version 0.0.1 2009.09.28 first version. """ from math import sqrt def gramm(X,inplace = False): # Returns the Gramm-Schmidt orthogonalization of matrix X if not inplace: V = [row[:] for row in X] # make a copy. else: V = X k = len(X[0]) # number of columns. n = len(X) # number of rows. for j in range(k): for i in range(j): # D = < Vi, Vj> D = sum([V[p][i]*V[p][j] for p in range(n)]) for p in range(n): # Note that the Vi's already have length one! # Vj = Vj - <Vi,Vj> Vi/< Vi,Vi > V[p][j] -= (D * V[p][i]) # Normalize column V[j] invnorm = 1.0 / sqrt(sum([(V[p][j])**2 for p in range(n)])) for p in range(n): V[p][j] *= invnorm return V |

The algorithm may be done in place or a separate matrix may be created. Note that the columns of the generated matrix V all have lengh or norm of one.

It is written so that it is stand-alone (does not require importing other Python modules). Here is the output for a random matrix with 10 columns.

Note: Have enclosed the code in CodeBox.

$ python gramm.py

i = 0 Input R:

0.71856 0.36612 0.69979 0.86844 0.25310 0.65541 0.11163 0.70527 0.63828 0.06808

0.88818 0.58432 0.75098 0.31272 0.61462 0.93716 0.60387 0.67479 0.21036 0.46674

0.69587 0.45944 0.54502 0.07974 0.61798 0.45362 0.13361 0.58535 0.93524 0.87348

0.95435 0.45627 0.88049 0.89684 0.35621 0.77600 0.53033 0.68745 0.28664 0.93202

0.87178 0.67467 0.83844 0.49449 0.88908 0.41148 0.70790 0.75325 0.43857 0.31378

0.28424 0.32494 0.22467 0.43489 0.97072 0.58225 0.67109 0.89471 0.52515 0.12127

0.80287 0.67120 0.74286 0.56414 0.16685 0.78518 0.75869 0.79553 0.16340 0.94957

0.24490 0.48795 0.66880 0.45138 0.84915 0.46690 0.08500 0.37051 0.36616 0.72665

0.82986 0.29793 0.81290 0.13830 0.30484 0.29833 0.83128 0.71961 0.64450 0.54384

0.10704 0.26563 0.00574 0.24027 0.61208 0.01169 0.88495 0.36639 0.15314 0.64091

Grammian G

G=

0.32311 -0.16905 0.14731 0.55100 -0.05067 -0.01341 -0.29920 0.31180 0.17256 -0.56786

0.39938 0.03371 -0.19987 -0.25739 0.00377 0.55781 0.05634 -0.49115 -0.06360 -0.41659

0.31291 0.02941 -0.25940 -0.38122 0.15226 0.02814 -0.57283 0.23063 0.45065 0.28323

0.42913 -0.27977 0.10755 0.45722 0.01883 -0.00606 0.01159 -0.44337 0.05061 0.56329

0.39201 0.21933 -0.04794 -0.13005 0.11465 -0.56237 -0.22895 -0.08496 -0.61773 -0.09155

0.12781 0.26481 -0.22145 0.27779 0.57380 0.37831 0.21333 0.40915 -0.26570 0.17736

0.36102 0.29383 -0.15635 0.00970 -0.73136 0.06587 0.28234 0.32969 -0.02801 0.17425

0.11012 0.61119 0.73796 -0.05701 0.09646 0.09787 -0.04031 -0.08804 0.19198 0.03715

0.37316 -0.42527 0.29062 -0.38205 0.22847 -0.17830 0.52688 0.25827 0.12790 -0.06949

0.04813 0.36356 -0.39329 0.17997 0.18732 -0.42750 0.34953 -0.23241 0.50433 -0.17338

D=

1.00000 -0.00000 -0.00000 -0.00000 0.00000 0.00000 -0.00000 0.00000 -0.00000 0.00000

-0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000 0.00000

-0.00000 0.00000 1.00000 -0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000 -0.00000

-0.00000 0.00000 -0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000

0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 -0.00000 -0.00000 -0.00000 0.00000

0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 -0.00000 0.00000

-0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000 1.00000 -0.00000 0.00000 -0.00000

0.00000 -0.00000 -0.00000 0.00000 -0.00000 0.00000 -0.00000 1.00000 0.00000 -0.00000

-0.00000 0.00000 0.00000 -0.00000 -0.00000 -0.00000 0.00000 0.00000 1.00000 0.00000

0.00000 0.00000 -0.00000 0.00000 0.00000 0.00000 -0.00000 -0.00000 0.00000 1.00000

$

Sep 28, 2009 | This entry needs more explanation. |

May 8, 1011 | The codebox above has been modified above to enable you to view the source code as saving to a file does not work anymore |

October 20th, 2010 at 12:47 pm

Thank you so much! I'm using your function here to help me display superhelices of an arbitrary nesting level in OpenGL

May 9th, 2011 at 11:13 pm

[...] The gramm.py code can be copied from from /p=184 [...]