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

3 Responses to “Gramm Schmidt orthogonalization algorithm in Python”

  1. Nathan Nifong Says:

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

  2. Solving least squares problem with the QR decomposition. · Digital explorations Says:

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

  3. Jordi Juanhuix Says:

    Thanks so much, I am using it to orthormalized measured reference systems of motorized instruments.

Leave a Reply