[Course] Advanced Computer Programming, Homework, week 12, NumPy

Some NumPy Exercises

Generate matrices A, with random Gaussian entries, B, a Toeplitz matrix, where $A∈R^{n×m}$ and $B∈R^{m×m}$, for n=200, m=500.

import numpy as np
from scipy.linalg import toeplitz

n, m = 200, 500
A = np.mat(np.random.randn(n,m))
A.shape

(200, 500)

B = np.mat(toeplitz(range(1,m+1), range(m+1,2*m+1)))
B.shape

(500, 500)

9.1 Matrix operations

A+A

matrix([[ 1.16129658e+00,  2.02802188e+00, -2.24653163e+00, ...,
          4.34751856e-01, -1.72832450e+00, -1.03742934e-01],
        [ 2.43605141e+00,  3.38559343e+00,  3.95117526e+00, ...,
          3.01734220e+00, -1.79109826e+00,  8.70547912e-01],
        [ 1.37962657e+00,  3.94905577e-01, -2.62851341e+00, ...,
         -2.30946756e+00,  6.35964393e-01, -1.77205804e+00],
        ...,
        [-2.26049647e+00, -1.88079505e+00, -3.52802465e-01, ...,
         -3.66617316e-01,  4.58773063e-01,  1.74914586e+00],
        [-3.17652236e-03, -7.45927965e-01, -5.21408396e-01, ...,
         -7.10983497e+00,  1.06145124e+00, -2.64221222e+00],
        [ 1.50903124e-01,  3.64857306e+00, -1.07906248e+00, ...,
          8.30657123e-01,  2.18301750e+00,  3.17541796e+00]])

A*A.T

matrix([[ 4.67285639e+02,  4.93627515e+00,  3.19671553e+01, ...,
          1.33838128e+00, -2.29387005e+01, -3.76985600e+01],
        [ 4.93627515e+00,  5.54675222e+02, -3.60075001e+01, ...,
         -4.62838913e+01,  3.22014590e+01, -5.12702264e-01],
        [ 3.19671553e+01, -3.60075001e+01,  5.28031417e+02, ...,
          1.16145442e+01,  2.44371450e+00, -5.87165498e+00],
        ...,
        [ 1.33838128e+00, -4.62838913e+01,  1.16145442e+01, ...,
          4.97031114e+02,  1.27083698e+01,  1.84652911e+01],
        [-2.29387005e+01,  3.22014590e+01,  2.44371450e+00, ...,
          1.27083698e+01,  5.33845828e+02,  1.63364502e+01],
        [-3.76985600e+01, -5.12702264e-01, -5.87165498e+00, ...,
          1.84652911e+01,  1.63364502e+01,  4.38938776e+02]])

A.T*A

matrix([[223.99945761,  -1.92497547,  19.25556243, ...,  15.22636245,
          16.78094612,   6.19667669],
        [ -1.92497547, 219.1314169 , -16.03395672, ...,   7.54911097,
         -12.82248639,  -2.24938105],
        [ 19.25556243, -16.03395672, 229.76560542, ..., -19.16528532,
          19.4542341 ,  16.23546023],
        ...,
        [ 15.22636245,   7.54911097, -19.16528532, ..., 206.23668558,
           2.59618655, -18.28650451],
        [ 16.78094612, -12.82248639,  19.4542341 , ...,   2.59618655,
         172.93390201,  -2.66838377],
        [  6.19667669,  -2.24938105,  16.23546023, ..., -18.28650451,
          -2.66838377, 181.05476695]])

A*B

matrix([[ -6768.88952192,  -6478.45424765,  -5969.30962591, ...,
           8165.39909155,   8276.96928964,   7846.04207362],
        [  1602.73655867,   2198.35168682,   3034.73791361, ...,
          13690.95765363,  14462.04753056,  14029.23619399],
        [  3937.1127593 ,   4269.55681991,   4356.21553895, ...,
          10749.73449171,  10187.34590325,  10361.95126696],
        ...,
        [    35.70041726,   -526.03502274,   -994.72590226, ...,
          -6135.59340453,  -6235.10482631,  -6127.80988023],
        [ -1362.21404422,  -1348.76015124,  -1521.74004678, ...,
         -10706.30883848, -12496.43802007, -12242.68419916],
        [ -1789.77364543,  -1750.66992742,   -833.50015227, ...,
           -974.77248501,   -773.69367354,   -232.34174922]])

def f(lamb):
    return A*(B-lamb*np.mat(np.eye(m,m)))

f(10)

matrix([[ -6774.69600481,  -6488.59435704,  -5958.07696774, ...,
           8163.22533227,   8285.61091215,   7846.56078829],
        [  1590.55630161,   2181.42371966,   3014.9820373 , ...,
          13675.87094262,  14471.00302185,  14024.88345443],
        [  3930.21462648,   4267.58229203,   4369.35810599, ...,
          10761.2818295 ,  10184.16608128,  10370.81155718],
        ...,
        [    47.0028996 ,   -516.6310475 ,   -992.96188994, ...,
          -6133.76031794,  -6237.39869162,  -6136.55560954],
        [ -1362.19816161,  -1345.03051142,  -1519.1330048 , ...,
         -10670.75966363, -12501.74527625, -12229.47313806],
        [ -1790.52816105,  -1768.91279272,   -828.10483986, ...,
           -978.92577062,   -784.60876106,   -248.21883899]])

9.2 Solving a linear system

b = np.random.rand(m)

x = np.linalg.solve(B,b)

9.3 Norms

forbenius_norm = np.linalg.norm(A, 'fro')
infinity_norm = np.linalg.norm(B, np.inf)
singular_values = np.linalg.svd(B)[1]
largest_singular_value = np.max(singular_values)
smallest_singualr_value = np.min(singular_values)

9.4 Power iteration

import numpy as np  
from scipy.linalg import toeplitz  
from time import clock

def power_iteration(A, n):  
    start = clock()
    u = np.random.randint(0, 20, size=n)
    num_iter = 0
    next_ = 1  
    p = 0 
    while num_iters==0 or abs(next_ - p) > 0.00001:  
        p = next_
        v = np.dot(A, u)  
        next_ = np.max(abs(v))  
        u = v/next_  
        num_iter += 1  
    return u, next_, num_iter, clock() - start  

np.random.seed(4)

for n in [50,100,200]:
    Z = np.random.normal(0.0, 0.1, (n, n))
    eigenvector, rou, num_iters, t = power_iteration(Z, n)  
    print('the largest eigenvalue: ', rou)  
    print('eigenvector: ', eigenvector)  
    print('the number of iterations: ', num_iters)  
    print('time: ', t)  

the largest eigenvalue:  0.7322776411310102
eigenvector:  [-0.37640648  0.19524665 -0.12680354 -0.66877869  0.09107082  0.38101643
 -0.25591548  0.35563192 -0.25498888  0.21865984  0.30670927  0.06145539
  0.2712122  -0.77382504 -0.09814322  0.40862076 -0.13278167 -0.00414074
 -0.45383608  0.15386053 -0.34115042 -0.11064735  0.17357024 -0.39248841
  0.11730325 -0.02913181 -1.         -0.56901765 -0.02391373 -0.15713455
 -0.28252059 -0.10739865  0.18981729  0.40356771  0.05537723  0.28917896
 -0.34084568 -0.33295133 -0.79586571  0.08282753  0.04320327 -0.13627579
  0.06495122  0.30480493  0.36914123 -0.37348225  0.03852317 -0.15413713
  0.28297282  0.03406147]
the number of iterations:  989
time:  0.013442279379887623
the largest eigenvalue:  1.3033154897500663
eigenvector:  [-0.01788451  0.14431714  0.27670934 -0.30243378  0.09606525  1.
  0.13820101  0.06744969 -0.01040092  0.29297793 -0.14568333  0.06115232
 -0.30641461  0.26253969  0.55194142 -0.34261743  0.287199    0.16196204
  0.04566519  0.40388928  0.11783789  0.33424901  0.10463527 -0.17237059
  0.12161516  0.25098784 -0.05471935 -0.49160685 -0.07928062 -0.17566036
  0.03879625 -0.26879562  0.27843559 -0.20557738 -0.35272519 -0.60193963
 -0.01222908 -0.19486227 -0.02153869  0.01396165  0.62101912 -0.22899374
 -0.29973264  0.08714527 -0.05255048 -0.18843342  0.16972285 -0.02569188
  0.28411925  0.149879    0.31977985  0.03451333 -0.28449795 -0.39621302
 -0.15735349  0.03662776 -0.13994955  0.30624533  0.0998112  -0.00602248
  0.26515036  0.05653403 -0.39760297  0.09622185  0.12490718  0.26016728
 -0.19233444  0.1310791  -0.02201675 -0.55281161 -0.15419681 -0.38813331
  0.41932544  0.02488147  0.23034242 -0.22447859  0.0227194  -0.54772363
 -0.15564906 -0.34916432 -0.19582291  0.32723005  0.28242736  0.25395502
  0.13588133  0.12638367  0.57853383 -0.18854702  0.11095108  0.20202662
  0.06579393 -0.10354477  0.29327317 -0.26299466  0.26865158 -0.63417166
  0.10887433 -0.20486772  0.36498914  0.21146259]
the number of iterations:  8807
time:  0.5758391682384172
the largest eigenvalue:  1.467981082770976
eigenvector:  [ 0.0381547   0.14024118 -0.47162136 -0.06915727  0.27727211 -0.49221961
 -0.43503035 -0.24934272  0.03467571 -0.19670185 -0.01232741  0.32816358
  0.30031234  0.14513348 -0.16837323  0.4151528  -0.11593887 -0.43600619
 -0.70835421  0.51850625  0.0315831  -0.31546986 -0.05815004 -0.48232464
  0.65186705 -0.05244794 -0.35929663 -0.13742196  0.41187342 -0.42661129
  0.64287751  0.23411644 -0.42253704 -0.07698448  0.14446508  0.22063569
  0.13292914  0.10264167 -0.20896717  0.43865298  0.08254183  0.01086428
  0.07315625 -0.10096209 -0.16480841 -0.36003873  0.74332686  0.1443572
 -0.36204752 -0.05779842 -0.27332622 -0.73065897  0.18959578  0.03113685
 -0.54223008 -0.24609017 -0.33823072  0.23486632 -0.08426762  0.25832281
 -0.54930581  0.08416751 -0.27978582  0.21354354 -0.46541108 -0.13185108
  0.07182352  0.51902111  0.44524258 -0.43795506  0.30087778 -0.24330084
  0.13242704  0.32727195 -0.32500576  0.10967076 -0.57448233  0.62779956
 -0.13778015  0.08650625 -0.02740524  0.66428279  0.37873053 -0.2069652
  0.55208693 -0.29206809 -0.17361356  0.71106908 -0.1987724   0.28627019
 -0.25716439  0.80297644 -0.67867617 -0.32622571  0.16976067 -0.04802861
  0.06264128  0.27858077 -0.31423294 -0.3997968  -0.2939908  -0.3065027
  0.15547628 -0.25683843  0.29837783 -0.06941952 -0.16303014 -0.31065087
 -0.58496087  0.21087573  0.08741621 -0.66290471 -0.43529586  0.04670688
 -0.4809971  -0.38116063 -0.32822356 -0.18952372 -0.33409959  0.08434364
  0.16314788 -0.4660477  -0.58817181  0.82948627 -0.00612028 -0.05741236
  0.51140909  0.00377282  0.04191805 -0.76875644  0.13814417  0.42758255
  0.36067425  0.22592804  0.194095    0.16453756 -0.2847642   0.5449104
  0.57491633  0.09508011  0.13828179 -0.02696916  0.36009175 -0.01068797
  0.25501144 -0.29714837 -0.52775785 -0.55497639 -0.21935517  0.29805999
  0.96601297 -0.25417463  0.31869312 -0.52237458 -0.09314444 -0.44935961
 -0.58632853  0.67977128  0.4525619   0.1951937  -1.          0.56672298
  0.14534756  0.27627565 -0.15444364 -0.43195298 -0.58118261 -0.5500504
 -0.39132198 -0.48566668  0.60482357 -0.15428003 -0.36395757  0.19977111
 -0.49051792  0.212237    0.02668716  0.19638928 -0.61574924  0.43274362
 -0.14769835 -0.09750915  0.28738431 -0.20951856  0.38691013 -0.24174646
 -0.05252725  0.09608743 -0.76198712 -0.25579623 -0.47036923  0.54259672
 -0.22032076  0.10842341 -0.16806426 -0.10727924 -0.13994377  0.1745847
  0.27877043  0.16804849]
the number of iterations:  604
time:  0.040909949414526636

9.5 Singular values

import scipy

for i in range(1, 200, 10):
    p = numpy.random.random()
    C = numpy.random.rand(i, i) > p
    rou = scipy.linalg.svd(C)[1].max()
    print(i, rou/(i*(1-p)))

1 0.0
11 1.368485097270882
21 1.4162962945755226
31 1.078834554574907
41 1.0102102029343
51 1.0549447195089567
61 1.0238686140130482
71 1.0146906824022468
81 0.9974517509985577
91 1.1323980627814303
101 1.0219351522124953
111 1.011452383838501
121 0.9990164076613034
131 1.002703317486608
141 1.0167804364212263
151 1.0960578967426344
161 1.0221028526150826
171 1.014370133299158
181 0.9976620627158513
191 1.002297239000159

可见，当n较大时，最大奇异值趋近于 $n(1-p)$

9.6 Nearest neighbor

def get_closest_v(A, z):
    return A[np.argmin(np.abs(A-z))]

A = np.random.rand(100)
get_closest_v(A, 0.5)

0.4981976220793418