orthogonal_matrices.py

orthogonal_matrices.py

>>> import numpy as np

>>> # No scientific notation
>>> np.set_printoptions(suppress=True)

>>> # Some random orthogonal matrix.
>>> orthog = np.array([[-0.3639,  0.8268,  0.3570,  0.2377], 
                       [-0.3578, -0.5330,  0.7427,  0.1906], 
                       [-0.6208, -0.1798, -0.5609,  0.5174], 
                       [-0.5951,  0.0024, -0.0797, -0.7997]])
>>> # Column one.
>>> orthog.T[0]
array([-0.3639, -0.3578, -0.6208, -0.5951])

>>> # Column two.
>>> orthog.T[1]
array([ 0.8268, -0.533 , -0.1798,  0.0024])

>>> # Lets assert the columns have a unit norm (excusing the
>>> # rounding error.)
>>> np.linalg.norm(orthog.T[0], ord=2)
0.99999034995343827
>>> np.linalg.norm(orthog.T[1], ord=2)
1.0000105199446654
>>> np.linalg.norm(orthog.T[2], ord=2)
1.0000065949782531
>>> np.linalg.norm(orthog.T[3], ord=2)
1.0000262496554777

>>> # Lets now assert that all the rows have a unit norm!
>>> np.linalg.norm(orthog[3], ord=2)
1.0000109749397752
>>> np.linalg.norm(orthog[2], ord=2)
1.0000161248699941
>>> np.linalg.norm(orthog[1], ord=2)
1.0000207447848271
>>> np.linalg.norm(orthog[0], ord=2)
0.99998586990017013

>>> # We can see that the column vectors or orthogonal to
>>> # each other as they return zero when we multiply them
>>> # together. This means the 'orthog' matrix is comprised
>>> # of orthonormal vectors.
>>> o = orthog.T[0].dot(orthog.T[1])
>>> '{:f}'.format(o)
'0.000026'

>>> # Print the results (slightly more) sensibly.
>>> np.set_printoptions(precision=1)

>>> # Here we get an identiy matrix (with some float rounding 
>>> # error!)
>>> orthog.dot(orthog.T)
array([[ 1., -0., -0., -0.],
       [-0.,  1., -0.,  0.],
       [-0., -0.,  1., -0.],
       [-0.,  0., -0.,  1.]])