In [84]:

    

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(42)


    

In [104]:

    

D = 16 # number of dimensions
condition = 1.e20 # approximate condition number we are going to make
scales = np.sqrt(condition) ** (np.arange(D) / D)
scales


    

    
Out[104]:





array([  1.00000000e+00,   4.21696503e+00,   1.77827941e+01,
         7.49894209e+01,   3.16227766e+02,   1.33352143e+03,
         5.62341325e+03,   2.37137371e+04,   1.00000000e+05,
         4.21696503e+05,   1.77827941e+06,   7.49894209e+06,
         3.16227766e+07,   1.33352143e+08,   5.62341325e+08,
         2.37137371e+09])

In [108]:

    

vectors = scales[:, None] * np.random.normal(size=(D, D+2))
matrix = np.dot(vectors, vectors.T)
matrix.shape


    

    
Out[108]:





(16, 16)

In [109]:

    

# this matrix has to be positive definite, but at high dynamic range, sometimes crazy shit happens!
eigvals, eigvecs = np.linalg.eigh(matrix)
"{:e}".format(np.max(eigvals) / np.min(np.abs(eigvals)))


    

    
Out[109]:





'2.001069e+19'

In [110]:

    

# hogg matrix inversion trick
foo = np.eye(D)
for iter in range(64): # no stopping criterion
    foo = np.dot(np.linalg.inv(np.dot(foo, matrix)), foo.T)
    foo = 0.5 * (foo + foo.T)
inverse = foo
np.allclose(np.dot(inverse, matrix), np.eye(D))


    

    
Out[110]:





False

In [ ]:

    




    

In [ ]: