from SciPy, scipy.linalg.svd - Singular Value Decomposition using SciPy

cf. scip.linalg.svd

Factorizes the matrix a into 2 unitary matrices U and Vh, and a 1-d. array s of singular values (real, non-negative) such that a == U*S*Vh, where S is a suitably shaped matrix of zeros with main diagonal s.


In [1]:
import numpy as np

In [2]:
from scipy import linalg

In [4]:
# Create an array of the given shape and populate it with
#    random samples from a uniform distribution
#    over ``[0, 1)``.
a = np.random.randn(9,6) + 1.j * np.random.randn(9,6)

In [5]:
a


Out[5]:
array([[ 0.32686765-1.03794508j,  0.65455117-0.88620527j,
        -0.22269114+0.34907433j, -0.85119968+0.3043245j ,
         0.35936646+0.01332564j,  1.57743118-1.62699021j],
       [-0.90979540-0.02522402j,  0.03367926-0.69622011j,
        -1.12272711+0.23931712j,  0.93699196-0.73548942j,
        -2.09994690+0.4551122j , -0.56650832-1.22022061j],
       [-0.17888326+1.51898273j,  1.48300148+0.71755258j,
         0.04284129+0.48819776j, -1.29268029+0.34231662j,
        -1.09280757-0.2872327j ,  1.00239153-0.11101026j],
       [-0.63645090-0.74049183j,  0.32565290+0.72510741j,
        -0.23952313+0.30741175j,  0.50011971+0.44605648j,
         1.34621709+0.12316658j, -1.36760350+2.46611406j],
       [-0.03144900+0.13232677j, -1.10060981+0.12429469j,
        -1.67518004+0.61615207j, -1.69533664-0.23601968j,
         1.14580207+0.30423478j,  1.57989044+0.47261328j],
       [-0.85501200-0.17574627j, -2.34780605+0.73568609j,
         1.19192757+0.06186622j, -1.28707194+0.92803125j,
        -0.42682367+0.84057498j, -0.37885884-0.18597101j],
       [ 0.75721205-1.48974891j, -0.78629419-0.70923246j,
        -0.27676234-1.23782299j, -0.01259970+1.35444433j,
        -1.28471679-0.25415455j,  0.94121300+1.13677475j],
       [-0.81240299+0.81865567j,  0.46555510+0.25448512j,
        -2.03706017-0.51214283j, -1.59398285+1.10851833j,
         1.21070496+0.81807185j,  0.51009566+0.58656106j],
       [ 1.18374747-0.9324303j ,  0.77354315-1.80342894j,
        -0.01615860-1.06586965j,  1.11528182+0.59689017j,
         0.96784141+1.07555765j, -3.06251784+0.18676989j]])

In [6]:
U, s, Vh = linalg.svd(a)

In [7]:
U.shape, Vh.shape, s.shape


Out[7]:
((9, 9), (6, 6), (6,))

In [8]:
U


Out[8]:
array([[-0.06974295-0.26101076j,  0.03950955-0.19081965j,
        -0.34330167+0.20153099j,  0.05392712+0.03927609j,
         0.00514477+0.45494835j,  0.26923882+0.3071468j ,
         0.31800010+0.1277802j , -0.14354608-0.16780999j,
        -0.31670108-0.29285557j],
       [-0.17795406-0.01140508j, -0.25325695-0.34254692j,
         0.08550203-0.18105876j,  0.42060028-0.0799726j ,
         0.21857863-0.39785314j, -0.10635291+0.12410679j,
         0.23180006+0.45134731j, -0.05769997+0.21981052j,
        -0.00083006-0.14649997j],
       [ 0.03934357-0.34783262j,  0.03009473+0.17141915j,
        -0.06690165-0.3277425j , -0.08129275-0.0229548j ,
         0.24268793-0.35866114j,  0.10945860+0.25387676j,
        -0.04719723-0.41813966j, -0.47332928-0.11215802j,
         0.09840385-0.20200591j],
       [ 0.20042190+0.34636829j,  0.17191821+0.30924989j,
         0.16207102+0.12187705j,  0.08501745-0.01589118j,
         0.05468756-0.09418463j,  0.21156484+0.61821912j,
        -0.15594416+0.2602303j , -0.13747514+0.03462962j,
        -0.19769000+0.27574045j],
       [ 0.35247775-0.12425068j, -0.18093378+0.06405496j,
        -0.41740149+0.24793054j,  0.04138898-0.09244324j,
        -0.04091100+0.06514755j, -0.31840098-0.09195336j,
         0.09183514+0.10906212j, -0.45992870+0.2739135j ,
         0.17697513+0.35151786j],
       [ 0.00152880+0.01980784j, -0.50241219+0.12280816j,
         0.11712159+0.26685802j, -0.41814537+0.23238973j,
        -0.07995375-0.1745932j ,  0.27049891-0.09891252j,
        -0.06255985+0.35934425j, -0.16744747-0.23024875j,
         0.20211067-0.20869028j],
       [ 0.07822782+0.07175871j,  0.12941850-0.30976238j,
        -0.14545936+0.1455782j , -0.19578789+0.62112746j,
        -0.06325809-0.38392166j, -0.12567985+0.17828758j,
         0.31634271-0.19958063j,  0.14641158-0.03001207j,
        -0.07159072+0.20784322j],
       [ 0.22346769-0.05730881j, -0.06569887+0.36693094j,
        -0.46141401+0.07730028j,  0.30115348+0.166406j  ,
         0.40146090-0.15521475j,  0.17617830-0.15017856j,
        -0.10691727+0.05050548j,  0.45118884-0.05194957j,
         0.01207156-0.10633819j],
       [-0.38390963+0.51047538j,  0.25320931+0.07040154j,
        -0.24382157+0.08438898j,  0.12637480+0.00604871j,
         0.07107945+0.05504163j,  0.09523111-0.00612539j,
         0.21389398-0.05465913j, -0.19070276-0.10138454j,
         0.56530613-0.09074874j]])

In [9]:
Vh


Out[9]:
array([[-0.26920770+0.j        , -0.20331352+0.19562025j,
        -0.26587793+0.02088816j, -0.16824448-0.21930761j,
         0.21158160-0.17991547j,  0.47558708+0.63013864j],
       [ 0.39385449+0.j        ,  0.52879684-0.2030625j ,
         0.04710054-0.05184119j,  0.33231677+0.22578786j,
         0.23571551-0.44854416j,  0.11675066+0.2954881j ],
       [-0.33936758+0.j        , -0.13773619+0.4759343j ,
         0.36178940+0.12456123j,  0.52669985-0.02691406j,
        -0.16632881-0.40098972j, -0.16217941+0.00636964j],
       [-0.31358909+0.j        ,  0.26832206+0.05588844j,
        -0.66869106+0.0141602j ,  0.50896971-0.01076258j,
         0.11081578+0.24253117j,  0.04440364-0.21668917j],
       [-0.40301820+0.j        ,  0.19108693-0.12905306j,
        -0.27909006-0.20178277j, -0.42068008+0.08104815j,
        -0.08306247-0.56178908j, -0.38492516-0.10763652j],
       [-0.62966158+0.j        ,  0.23598612-0.41239851j,
         0.45980278+0.01360805j,  0.01169945+0.20298486j,
         0.14460131+0.25126593j,  0.18136226+0.08879425j]])

In [10]:
s


Out[10]:
array([ 6.49313842,  4.91573534,  4.21719281,  3.8309913 ,  2.48135155,
        1.72791536])

full_matrices : bool, optional - if True, U and Vh are of shape (M,M), (N,N). If False, the shapes are (M,K) and (K,N) where K = min(M,N)


In [11]:
U,s,Vh=linalg.svd(a,full_matrices=False)

In [12]:
U.shape, Vh.shape, s.shape


Out[12]:
((9, 6), (6, 6), (6,))

In [13]:
S=linalg.diagsvd(s,6,6)  # Construct the Sigma matrix S, given the vector s, of a == U*S*Vh

In [14]:
S


Out[14]:
array([[ 6.49313842,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  4.91573534,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  4.21719281,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  3.8309913 ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  2.48135155,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         1.72791536]])

In [15]:
np.allclose(a,np.dot(U,np.dot(S,Vh)))


Out[15]:
True

compute_uv : bool, optional Whether to compute also U and Vh in addition to s. Default is True.


In [16]:
s2 = linalg.svd(a,compute_uv=False)

In [17]:
np.allclose(s,s2)


Out[17]:
True

More Simple examples

"Gold standard" example for cuSOLVER in CUDA Toolkit Documentation


In [3]:
A_gold = np.array( [[1.0,2.0], [4.0,5.0],[2.,1.0]])

In [5]:
UsVh_gold = linalg.svd(A_gold)

In [6]:
UsVh_gold[0]


Out[6]:
array([[-0.30821892,  0.48819507,  0.81649658],
       [-0.90613333,  0.11070553, -0.40824829],
       [-0.28969549, -0.86568462,  0.40824829]])

In [7]:
linalg.diagsvd(UsVh_gold[1],2,2)


Out[7]:
array([[ 7.0652835,  0.       ],
       [ 0.       ,  1.0400813]])

In [8]:
UsVh_gold[2]


Out[8]:
array([[-0.63863584, -0.76950911],
       [-0.76950911,  0.63863584]])

In [5]:
A_6_241_04 = np.array([[100,100],[100.2,100]])
A_6_241_04


Out[5]:
array([[ 100. ,  100. ],
       [ 100.2,  100. ]])

In [7]:
UsVh_6_241_04 = linalg.svd(A_6_241_04)

In [8]:
UsVh_6_241_04[0]


Out[8]:
array([[-0.70675314, -0.70746025],
       [-0.70746025,  0.70675314]])

In [9]:
linalg.diagsvd( UsVh_6_241_04[1], 2,2)


Out[9]:
array([[  2.00100050e+02,   0.00000000e+00],
       [  0.00000000e+00,   9.99500000e-02]])

In [10]:
UsVh_6_241_04[2]


Out[10]:
array([[-0.70746025, -0.70675314],
       [ 0.70675314, -0.70746025]])

In [11]:
A33 = np.array([[-149 ,-50, -154 ],[537 ,180, 546 ],[-27 ,-9, -25  ]])

In [13]:
UsVh_33 = linalg.svd(A33)

In [14]:
UsVh_33[0]


Out[14]:
array([[-0.26906707, -0.67982121,  0.68223605],
       [ 0.96200923, -0.15566953,  0.22428829],
       [-0.04627257,  0.71666597,  0.69587983]])

In [15]:
linalg.diagsvd( UsVh_33[1], 3,3)


Out[15]:
array([[  8.17759668e+02,   0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   2.47497449e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   2.96452308e-03]])

In [16]:
UsVh_33[2]


Out[16]:
array([[ 0.68227785,  0.22871202,  0.6943974 ],
       [-0.66714135, -0.19371852,  0.71930213],
       [ 0.29903068, -0.95402513,  0.02041339]])

Simple Examples, Test cases, of applying SVD successfully to obtain a Matrix Product State (MPS)

Let $d=2$ (dimension of the state space, for each, say, spin system, $L=2$ (number of sites; should be an even number)


In [52]:
d=2
L=2
Psi_d2_L2 = np.random.randn(d**(L-1),d) + 1.j * np.random.randn(d**(L-1),d)

In [53]:
print(Psi_d2_L2.shape)
Psi_d2_L2


(2, 2)
Out[53]:
array([[  1.63401174e+00-0.13735011j,  -7.01045042e-01+2.08744185j],
       [ -1.27862833e-03+0.42564405j,  -9.93358872e-01+0.83671748j]])

In [54]:
UsVh_d2_L2 = linalg.svd(Psi_d2_L2)

In [55]:
print(UsVh_d2_L2[0].shape)
UsVh_d2_L2[0]


(2, 2)
Out[55]:
array([[-0.91091206-0.0026881j ,  0.31731043-0.26371592j],
       [-0.30341964-0.27958633j, -0.90985930-0.04386417j]])

In [56]:
UsVh_d2_L2[1]


Out[56]:
array([ 3.00018954,  0.63610712])

In [57]:
linalg.diagsvd(UsVh_d2_L2[1],2,2)


Out[57]:
array([[ 3.00018954,  0.        ],
       [ 0.        ,  0.63610712]])

In [58]:
print(UsVh_d2_L2[2].shape)
UsVh_d2_L2[2]


(2, 2)
Out[58]:
array([[-0.53552886+0.j        ,  0.23346824-0.81160423j],
       [ 0.84451693+0.j        ,  0.14804792-0.51465811j]])

Calculate US and reshape


In [59]:
US_d2_L2 = np.dot( UsVh_d2_L2[0] , linalg.diagsvd(UsVh_d2_L2[1],2,2) )
print(US_d2_L2.shape)
US_d2_L2


(2, 2)
Out[59]:
array([[-2.73290882-0.00806481j,  0.20184343-0.16775157j],
       [-0.91031644-0.83881198j, -0.57876798-0.02790231j]])

In [63]:
# new matrix after 1 iteration, l=1
Psi_d2_L2_l1 = US_d2_L2.reshape(d**(L-(1+1)),d*2)
print(Psi_d2_L2_l1.shape)
Psi_d2_L2_l1


(1, 4)
Out[63]:
array([[-2.73290882-0.00806481j,  0.20184343-0.16775157j,
        -0.91031644-0.83881198j, -0.57876798-0.02790231j]])

Calculate the right normalized matrices $B$'s which are only columns of Vh in this case:


In [64]:
B0_d2_L2=UsVh_d2_L2[2][:,0]
B1_d2_L2=UsVh_d2_L2[2][:,1]
print(B0_d2_L2)
print(B1_d2_L2)


[-0.53552886+0.j  0.84451693+0.j]
[ 0.23346824-0.81160423j  0.14804792-0.51465811j]

Examples with complex numbers


In [33]:
M=4
N=2
ind_RR = 1
ind_CC = 0.1
A_CC=[]
for row in range(M):
    A_row =[]
    for col in range(N):
        A_val = ind_RR * (row+1 + M*col) + ind_CC*(row+1 + M*col)*1j
        A_row.append(A_val)
    A_CC.append(A_row)

In [34]:
A_CC = np.array(A_CC)
A_CC


Out[34]:
array([[ 1.+0.1j,  5.+0.5j],
       [ 2.+0.2j,  6.+0.6j],
       [ 3.+0.3j,  7.+0.7j],
       [ 4.+0.4j,  8.+0.8j]])

In [35]:
UsVh_CC = linalg.svd(A_CC)

In [36]:
print(UsVh_CC[0].shape)
UsVh_CC[0]


(4, 4)
Out[36]:
array([[-0.35031448-0.03503145j,  0.75521459+0.07552146j,
        -0.39607577+0.04439239j, -0.36708362+0.08001822j],
       [-0.44142415-0.04414242j,  0.31964734+0.03196473j,
         0.24825442-0.06940632j,  0.78585375-0.1264373j ],
       [-0.53253383-0.05325338j, -0.11591991-0.01159199j,
         0.69171847+0.00563549j, -0.47045664+0.01281994j],
       [-0.62364350-0.06236435j, -0.55148716-0.05514872j,
        -0.54389712+0.01937845j,  0.05168651+0.03359914j]])

In [37]:
print(UsVh_CC[1].shape)
linalg.diagsvd(UsVh_CC[1],2,2)


(2,)
Out[37]:
array([[ 14.29836749,   0.        ],
       [  0.        ,   1.26360085]])

In [38]:
print(UsVh_CC[2].shape)
UsVh_CC[2]


(2, 2)
Out[38]:
array([[-0.37616823 +0.00000000e+00j, -0.92655138 +4.00453083e-18j],
       [-0.92655138 +0.00000000e+00j,  0.37616823 -1.62578927e-18j]])

$d=2$, $L=4$


In [82]:
d=2
L=4
Psi_d2_L4 = np.random.randn(d**(L-1),d) + 1.j * np.random.randn(d**(L-1),d)

In [83]:
print(Psi_d2_L4.shape)
Psi_d2_L4


(8, 2)
Out[83]:
array([[-0.54089591+0.24894326j, -0.77157693+0.59605689j],
       [-1.28483134-0.88204462j, -1.57853063-0.03664469j],
       [-0.67814438-0.52990602j,  0.83399876-1.13522893j],
       [-0.08278651-0.84136502j,  0.27905342-0.51229099j],
       [ 0.38046230-0.2488228j , -0.20882187-0.55554895j],
       [-0.32651707-0.14839783j, -0.73316458-0.75729396j],
       [ 0.82195234+0.20466075j, -0.21694321-1.06020886j],
       [ 2.88306483+0.45530495j,  0.36009038+0.76975469j]])

In [84]:
UsVh_d2_L4 = linalg.svd(Psi_d2_L4)

In [85]:
print(UsVh_d2_L4[0].shape)
UsVh_d2_L4[0]


(8, 8)
Out[85]:
array([[ -1.90485504e-01 +1.50115593e-01j,
          9.25232955e-02 -2.08046244e-01j,
         -2.16395809e-01 +1.64131486e-01j,
         -7.95721023e-03 +2.31934285e-01j,
          1.10238758e-01 +3.48780919e-02j,
         -8.26533455e-02 +1.83355013e-02j,
          2.04870519e-01 -1.24073021e-01j,
          7.98215402e-01 -2.16016105e-01j],
       [ -4.65948857e-01 -1.39844728e-01j,
          2.55479921e-01 -3.01641527e-01j,
          3.77495107e-01 +2.53405485e-01j,
          7.97451425e-02 +8.71903084e-02j,
         -8.39276821e-02 +2.09721958e-01j,
         -1.89595118e-01 +2.73723823e-01j,
         -6.07701593e-02 +4.49126676e-01j,
         -1.22709419e-01 -6.35263800e-02j],
       [ -7.32053830e-02 -2.77826445e-01j,
         -3.22827451e-01 +3.32325361e-01j,
          7.29332352e-01 +1.47531924e-02j,
         -8.20093368e-02 +5.59982886e-02j,
          1.00629398e-02 -6.24653596e-02j,
         -2.34090026e-02 -1.40098276e-01j,
          1.18085815e-03 -1.71434953e-01j,
          3.24146818e-01 +4.63666337e-02j],
       [  1.50626381e-03 -2.48187634e-01j,
         -6.70809833e-02 +1.64352483e-02j,
         -8.71245595e-02 -5.45799386e-02j,
          9.40352261e-01 +1.43079033e-03j,
          1.50357243e-04 +2.33366854e-03j,
          4.36928205e-03 -4.82849476e-02j,
          2.96973675e-02 -6.96122087e-03j,
          9.27779517e-02 +1.64019624e-01j],
       [  3.62714790e-02 -1.14610279e-01j,
          1.93833867e-01 +1.14059768e-01j,
          2.35463291e-02 +5.49029275e-02j,
         -2.05959637e-03 -3.45313006e-03j,
          9.56702824e-01 +5.33343457e-03j,
         -4.31870528e-02 +2.55475863e-04j,
         -7.95803807e-02 +4.42260322e-02j,
         -7.68715727e-02 -3.25481480e-03j],
       [ -1.82935415e-01 -1.02440610e-01j,
          2.35891426e-01 +1.58345591e-01j,
         -8.14403150e-03 +1.35067051e-01j,
          4.35688095e-03 +4.70321384e-02j,
         -4.33688298e-02 +9.26384233e-03j,
          9.12914834e-01 +5.08308989e-03j,
         -8.88893105e-02 +2.95481412e-02j,
          4.88048569e-02 -9.06625803e-02j],
       [  1.11839314e-01 -8.23779530e-02j,
          3.30052296e-01 +3.74244104e-01j,
          3.60284252e-02 +1.66936958e-01j,
          2.67068939e-02 +1.76183818e-03j,
         -8.88210907e-02 -2.31184747e-02j,
         -9.59071670e-02 -1.01023574e-02j,
          7.95117698e-01 +2.51739085e-02j,
         -1.62567573e-01 -1.57582764e-01j],
       [  6.71097852e-01 +1.79701690e-01j,
          4.18917513e-01 -1.29176530e-01j,
          3.45190684e-01 -8.12589702e-02j,
          8.09441485e-02 -1.68939289e-01j,
         -7.16288743e-02 +2.76652707e-02j,
          5.74116388e-02 +1.07397863e-01j,
         -1.27511521e-01 +2.01831598e-01j,
          2.91477963e-01 +4.02915065e-02j]])

In [86]:
UsVh_d2_L4[1]


Out[86]:
array([ 4.07925268,  2.50908488])

In [87]:
print(UsVh_d2_L4[2].shape)
UsVh_d2_L4[2]


(2, 2)
Out[87]:
array([[ 0.85234517+0.j        ,  0.50740978+0.12666106j],
       [ 0.52297965+0.j        , -0.82696961-0.20643049j]])

In [88]:
US_d2_L4 = np.dot( UsVh_d2_L4[0] , linalg.diagsvd(UsVh_d2_L4[1], UsVh_d2_L4[0].shape[0],d) )
print(US_d2_L4.shape)
US_d2_L4


(8, 2)
Out[88]:
array([[-0.77703850+0.61235944j,  0.23214880-0.52200569j],
       [-1.90072312-0.57046198j,  0.64102081-0.75684419j],
       [-0.29862325-1.13332427j, -0.81000148+0.83383254j],
       [ 0.00614443-1.01242007j, -0.16831188+0.04123743j],
       [ 0.14796053-0.46752429j,  0.48634562+0.28618564j],
       [-0.74623978-0.41788113j,  0.59187161+0.39730253j],
       [ 0.45622082-0.33604049j,  0.82812923+0.93901022j],
       [ 2.73757771+0.7330486j ,  1.05109960-0.32411488j]])

In [89]:
# after iteration l = 1, we obtain a new matrix \Psi to apply SVD on
l=1
Psi_d2_L4_l1 = US_d2_L4.reshape(d**(L-(l+1)), d*2 )
print(Psi_d2_L4_l1.shape )


(4, 4)

In [90]:
UsVh_d2_L4_l2 = linalg.svd(Psi_d2_L4_l1)

In [91]:
print(UsVh_d2_L4_l2[0].shape)
UsVh_d2_L4_l2[0]


(4, 4)
Out[91]:
array([[-0.50481896+0.22658928j,  0.28752163-0.32264486j,
         0.08050191-0.37960522j,  0.10490858+0.58775865j],
       [-0.12666083-0.22143248j, -0.09077695+0.62361639j,
        -0.45751654-0.55642761j,  0.13028581+0.04335185j],
       [-0.15186365-0.02584841j,  0.32851475+0.20651011j,
         0.52211188-0.09498952j,  0.61310101-0.41010273j],
       [ 0.64575814-0.43359604j,  0.44438420-0.26086024j,
        -0.09834857-0.19796462j,  0.18725103+0.21340134j]])

In [92]:
print(UsVh_d2_L4_l2[1].shape)
UsVh_d2_L4[1]


(4,)
Out[92]:
array([ 4.07925268,  2.50908488])

In [93]:
print(UsVh_d2_L4_l2[2].shape)
UsVh_d2_L4_l2[2]


(4, 4)
Out[93]:
array([[ 0.31233239+0.j        , -0.06776646+0.21505645j,
         0.65673164+0.63799897j,  0.05903944+0.09890462j],
       [-0.38260696+0.j        ,  0.52190038+0.43742359j,
        -0.13397096+0.16384349j,  0.57852121+0.10204319j],
       [ 0.51863095+0.j        ,  0.06982417-0.43719057j,
        -0.11813854-0.10114526j,  0.52662728+0.48320551j],
       [ 0.69791284+0.j        ,  0.26455367+0.46844356j,
        -0.27955709-0.12053514j, -0.10061270-0.34739869j]])

In [94]:
US_d2_L4_l2 = np.dot( UsVh_d2_L4_l2[0] , linalg.diagsvd(UsVh_d2_L4_l2[1], UsVh_d2_L4_l2[0].shape[1],UsVh_d2_L4_l2[2].shape[0]) )
print(US_d2_L4_l2.shape)
US_d2_L4_l2


(4, 4)
Out[94]:
array([[-2.01991593+0.90664445j,  0.64490984-0.72369109j,
         0.09554528-0.45054196j,  0.07313365+0.40973708j],
       [-0.50680391-0.88601068j, -0.20361232+1.398769j  ,
        -0.54301254-0.6604071j ,  0.09082457+0.03022135j],
       [-0.60764717-0.1034264j ,  0.73685725+0.46320133j,
         0.61967880-0.11274019j,  0.42740369-0.28588996j],
       [ 2.58385138-1.73493397j,  0.99675194-0.58510844j,
        -0.11672694-0.23495823j,  0.13053605+0.14876589j]])

In [97]:
# after iteration l = 2, we obtain a new matrix \Psi to apply SVD on
l=2
Psi_d2_L4_l2 = US_d2_L4_l2.reshape(d**(L-(l+1)), d* US_d2_L4_l2.shape[1] )
print(Psi_d2_L4_l2.shape )


(2, 8)

In [98]:
UsVh_d2_L4_l3 = linalg.svd(Psi_d2_L4_l2)

In [99]:
print(UsVh_d2_L4_l3[0].shape)
UsVh_d2_L4_l3[0]


(2, 2)
Out[99]:
array([[-0.56219326 -3.08515691e-19j,  0.82700589 -2.82510600e-18j],
       [-0.16549764 -8.10277279e-01j, -0.11250423 -5.50821260e-01j]])

In [96]:
L-(l+1)


Out[96]:
1

Examples with fixed values


In [60]:
d=2
L=2

In [3]:
def create_fixed_CC_mat(d,L):
    totalsysstates = d**(L-1)
    A = []
    for i in range(totalsysstates):
        ithstate = []
        f = i*(0.9/totalsysstates)+0.1
        theta_f = 2.*np.arccos(-1.)*f  
        d0 = f*(np.cos( theta_f) + np.sin(theta_f)*1j)
        d1 = (1.-f)*(np.sin( theta_f) + np.cos(theta_f)*1j)
        ithstate=[d0,d1]
        A.append(ithstate)
    return np.array(A)

In [66]:
A_CC_d2L2=create_fixed_CC_mat(d,L)
print(A_CC_d2L2.shape)
print(A_CC_d2L2)


(2, 2)
[[ 0.08090170+0.05877853j  0.52900673+0.72811529j]
 [-0.52308108-0.16995935j -0.13905765-0.42797543j]]

In [67]:
UsVh_CC_d2L2 = linalg.svd(A_CC_d2L2)

In [68]:
print(UsVh_CC_d2L2[0].shape)
print(UsVh_CC_d2L2[0])


(2, 2)
[[-0.80191235-0.14797714j -0.57881967+0.00267245j]
 [ 0.52737106+0.23857728j -0.79251382-0.19204816j]]

In [69]:
print(UsVh_CC_d2L2[1].shape)
print(linalg.diagsvd(UsVh_CC_d2L2[1],2,2 ) )


(2,)
[[ 1.06765472  0.        ]
 [ 0.          0.43024808]]

In [70]:
print(UsVh_CC_d2L2[2].shape)
print(UsVh_CC_d2L2[2])


(2, 2)
[[-0.36526811+0.j         -0.66257530-0.65389081j]
 [ 0.93090236+0.j         -0.25998175-0.25657412j]]

In [71]:
Psi_new_CC_d2L2 = np.dot( UsVh_CC_d2L2[0], linalg.diagsvd(UsVh_CC_d2L2[1],2,2 ) )
print(Psi_new_CC_d2L2)
Psi_new_CC_d2L2 = Psi_new_CC_d2L2.reshape(1,4)
print(Psi_new_CC_d2L2)


[[-0.85616550-0.15798849j -0.24903605+0.00114982j]
 [ 0.56305020+0.25471816j -0.34097755-0.08262835j]]
[[-0.85616550-0.15798849j -0.24903605+0.00114982j  0.56305020+0.25471816j
  -0.34097755-0.08262835j]]
$d=2,L=4$

In [4]:
d=2
L=4
A_CC_d2L4=create_fixed_CC_mat(d,L)
print(A_CC_d2L4.shape)
print(A_CC_d2L4)


(8, 2)
[[ 0.08090170+0.05877853j  0.52900673+0.72811529j]
 [ 0.04960714+0.20662861j  0.76574131+0.18383822j]
 [-0.14754691+0.28957712j  0.60142940-0.30644359j]
 [-0.40419730+0.167424j    0.21525943-0.51968224j]
 [-0.52308108-0.16995935j -0.13905765-0.42797543j]
 [-0.34615530-0.56487411j -0.28776606-0.17634327j]
 [ 0.12123671-0.76545846j -0.22222988+0.03519775j]
 [ 0.67486029-0.57638514j -0.07306291+0.08554567j]]

In [5]:
UsVh_CC_d2L4 = linalg.svd(A_CC_d2L4)

In [6]:
print(UsVh_CC_d2L4[0].shape)
print(UsVh_CC_d2L4[0])


(8, 8)
[[-0.33538340 +2.94252659e-01j  0.30493681 -1.66022303e-01j
   0.08299306 -3.37201981e-02j  0.30697852 -2.27889872e-02j
   0.38434096 -9.57592973e-02j  0.25136223 -2.32029589e-01j
  -0.04234271 -3.28344766e-01j -0.35084081 -2.66462847e-01j]
 [-0.06891648 +4.60033722e-01j  0.08846466 -1.69795370e-01j
  -0.44825252 -2.43956374e-01j -0.22278302 -2.45608964e-01j
   0.05147936 -6.19234440e-02j  0.27017747 +1.31082756e-01j
   0.37038487 +1.62678423e-01j  0.34069534 -1.32750330e-02j]
 [ 0.10941525 +4.01178272e-01j -0.19426785 -6.24947677e-02j
   0.82039936 +1.04838159e-02j -0.13751941 +6.36205675e-03j
  -0.02311415 +3.96304010e-02j  0.10854465 +6.09296957e-02j
   0.19115862 +2.32716587e-02j  0.18050863 -7.78194422e-02j]
 [ 0.12517450 +1.61900678e-01j -0.40641482 +1.43365508e-02j
  -0.13598571 +5.79701065e-02j  0.83937790 +3.55872198e-02j
  -0.09660613 +4.30692460e-02j  0.02246140 +5.89198402e-02j
   0.13240239 +4.37328096e-02j  0.17262781 -2.56815645e-02j]
 [ 0.03460551 -1.32791451e-01j -0.43416981 -3.72932078e-02j
  -0.04074486 +3.90318696e-02j -0.11652017 +6.26847294e-05j
   0.86075387 -4.84603034e-03j -0.09717149 +3.30587115e-02j
  -0.00867450 +8.26112532e-02j  0.08568918 +9.38919831e-02j]
 [-0.03288777 -3.41958367e-01j -0.24712181 -1.95553782e-01j
   0.05803209 -1.03053653e-02j -0.02008655 -6.84519260e-02j
  -0.11400378 -8.86210720e-02j  0.82376872 -3.42350796e-02j
  -0.16347369 +8.01811866e-02j -0.06383315 +1.90432380e-01j]
 [ 0.02203726 -3.76133912e-01j  0.07897272 -3.27859171e-01j
   0.11822228 -3.38449394e-02j  0.09040788 -1.06594482e-01j
  -0.01693170 -1.57281628e-01j -0.15716687 -1.28425828e-01j
   0.75174616 -2.66339805e-03j -0.21916242 +1.73985639e-01j]
 [ 0.18895214 -2.35038590e-01j  0.39184287 -2.81594907e-01j
   0.12135809 +2.25227404e-03j  0.16782898 -6.43926928e-02j
   0.11325531 -1.53989533e-01j -0.03600064 -2.01339328e-01j
  -0.20884118 -1.46372454e-01j  0.69670082 +1.91546343e-02j]]

In [7]:
print(UsVh_CC_d2L4[1].shape)
print(linalg.diagsvd(UsVh_CC_d2L4[1],2,2 ) )


(2,)
[[ 1.63372423  0.        ]
 [ 0.          1.54748025]]

In [8]:
print(UsVh_CC_d2L4[2].shape)
print(UsVh_CC_d2L4[2])


(2, 2)
[[ 0.56369984+0.j          0.01471577-0.82584862j]
 [ 0.82597972+0.j         -0.01004296+0.56361037j]]

In [18]:
Psi_new_CC_d2L4 = np.dot( UsVh_CC_d2L4[0], linalg.diagsvd(UsVh_CC_d2L4[1],8,2 ) )
print(Psi_new_CC_d2L4)
Psi_new_CC_d2L4 = Psi_new_CC_d2L4.reshape( (2**(L-2),d*d),order='F')
print(Psi_new_CC_d2L4)


[[-0.54792399+0.4807277j   0.47188370-0.25691624j]
 [-0.11259053+0.75156824j  0.13689731-0.26275498j]
 [ 0.17875434+0.65541466j -0.30062567-0.09670942j]
 [ 0.20450062+0.26450106j -0.62891891+0.02218553j]
 [ 0.05653586-0.21694461j -0.67186921-0.0577105j ]
 [-0.05372954-0.55866567j -0.38241613-0.30261562j]
 [ 0.03600281-0.61449909j  0.12220873-0.50735559j]
 [ 0.30869569-0.38398824j  0.60636911-0.43576256j]]
[[-0.54792399+0.4807277j   0.05653586-0.21694461j  0.47188370-0.25691624j
  -0.67186921-0.0577105j ]
 [-0.11259053+0.75156824j -0.05372954-0.55866567j  0.13689731-0.26275498j
  -0.38241613-0.30261562j]
 [ 0.17875434+0.65541466j  0.03600281-0.61449909j -0.30062567-0.09670942j
   0.12220873-0.50735559j]
 [ 0.20450062+0.26450106j  0.30869569-0.38398824j -0.62891891+0.02218553j
   0.60636911-0.43576256j]]

In [19]:
UsVh_CC_d2L4l01 = linalg.svd(Psi_new_CC_d2L4)

In [20]:
print(UsVh_CC_d2L4l01[0].shape)
print(UsVh_CC_d2L4l01[0])


(4, 4)
[[-0.18535573+0.33672088j -0.61731879+0.16726552j  0.50657445+0.34085335j
  -0.08184368+0.25239847j]
 [ 0.00467968+0.61092099j -0.25600887+0.05655513j -0.35481630-0.05715464j
   0.14777127-0.63797885j]
 [ 0.04464389+0.60308742j  0.23327541-0.14688632j -0.35310495-0.06081485j
  -0.14997130+0.63829972j]
 [-0.08528068+0.32563115j  0.63012779-0.22169054j  0.50331445+0.34107634j
   0.08579891-0.25192038j]]

In [21]:
print(UsVh_CC_d2L4l01[1].shape)
print(linalg.diagsvd(UsVh_CC_d2L4l01[1],4,4 ) )


(4,)
[[ 1.68907122  0.          0.          0.        ]
 [ 0.          1.48037668  0.          0.        ]
 [ 0.          0.          0.13718988  0.        ]
 [ 0.          0.          0.          0.02126598]]

In [22]:
print(UsVh_CC_d2L4l01[2].shape)
print(UsVh_CC_d2L4l01[2])


(4, 4)
[[ 0.70689517+0.j         -0.55973645-0.03879965j -0.20410253+0.10879021j
  -0.34083506+0.12580633j]
 [ 0.34155756+0.j          0.19540716-0.02773377j -0.56831742-0.03578054j
   0.72118017+0.00442023j]
 [-0.19351374+0.j          0.39906275-0.09767257j -0.57603008+0.37616081j
  -0.45754797+0.33332217j]
 [ 0.58837920+0.j          0.69029687+0.03059078j  0.38567384+0.01378384j
  -0.15964484-0.0440859j ]]

In [26]:
np.dot( UsVh_CC_d2L4l01[0], linalg.diagsvd(UsVh_CC_d2L4l01[1],4,4 ) )


Out[26]:
array([[-0.31307902+0.56874555j, -0.91386435+0.24761598j,
         0.06949689+0.04676163j, -0.00174049+0.0053675j ],
       [ 0.00790432+1.03188905j, -0.37898955+0.0837229j ,
        -0.04867720-0.00784104j,  0.00314250-0.01356725j],
       [ 0.07540670+1.0186576j ,  0.34533548-0.21744708j,
        -0.04844242-0.00834318j, -0.00318929+0.01357407j],
       [-0.14404515+0.5500142j ,  0.93282648-0.32818551j,
         0.06904965+0.04679222j,  0.00182460-0.00535733j]])

In [24]:
Psi_new_CC_d2L4l01 = np.dot( UsVh_CC_d2L4l01[0], linalg.diagsvd(UsVh_CC_d2L4l01[1],4,4 ) )
print(Psi_new_CC_d2L4l01)
Psi_new_CC_d2L4l01 = Psi_new_CC_d2L4l01.reshape( (2**(L-3),d*d*d),order='F')
print(Psi_new_CC_d2L4l01)


[[-0.31307902+0.56874555j -0.91386435+0.24761598j  0.06949689+0.04676163j
  -0.00174049+0.0053675j ]
 [ 0.00790432+1.03188905j -0.37898955+0.0837229j  -0.04867720-0.00784104j
   0.00314250-0.01356725j]
 [ 0.07540670+1.0186576j   0.34533548-0.21744708j -0.04844242-0.00834318j
  -0.00318929+0.01357407j]
 [-0.14404515+0.5500142j   0.93282648-0.32818551j  0.06904965+0.04679222j
   0.00182460-0.00535733j]]
[[-0.31307902+0.56874555j  0.07540670+1.0186576j  -0.91386435+0.24761598j
   0.34533548-0.21744708j  0.06949689+0.04676163j -0.04844242-0.00834318j
  -0.00174049+0.0053675j  -0.00318929+0.01357407j]
 [ 0.00790432+1.03188905j -0.14404515+0.5500142j  -0.37898955+0.0837229j
   0.93282648-0.32818551j -0.04867720-0.00784104j  0.06904965+0.04679222j
   0.00314250-0.01356725j  0.00182460-0.00535733j]]

Calculate the dimensions

$L=2$ case


In [100]:
d=2
L=2

In [101]:
d**L


Out[101]:
4

In [102]:
d**(L-1)


Out[102]:
2

In [103]:
d**(L/2)


Out[103]:
2

In [104]:
L/2


Out[104]:
1

In [107]:
[l+1 for l in range(L/2)]


Out[107]:
[1]

In [108]:
range(1,2)


Out[108]:
[1]

In [109]:
def calculate_dims(d,L):
    results = []
    result_1 = []
    result_1.append( (d**(L-1),d) )
    r_1 = min( result_1[0][0], result_1[0][1])
    newPsidim = ( result_1[0][0]/d, d*r_1) 
    result_1.append( newPsidim )
    newBdim   = ( r_1, 1)
    result_1.append( newBdim )
    results.append( result_1 )
    for l in range(2,L+1):
        result_l = []
        r_previous = results[l-2][2][0]
        result_l.append( (d**(L-l) , d*r_previous))
        r_l = min( result_l[0][0], result_l[0][1])
        newPsidim = ( result_l[0][0]/d , d*r_l)
        result_l.append( newPsidim )
        newBdim = ( r_l, r_previous)
        result_l.append( newBdim)
        results.append( result_l )
    return results

In [110]:
results_d2_L2 = calculate_dims(2,2)

In [111]:
results_d2_L2


Out[111]:
[[(2, 2), (1, 4), (2, 1)], [(1, 4), (0, 2), (1, 2)]]

In [112]:
results_d2_L4 = calculate_dims(2,4)

In [113]:
results_d2_L4


Out[113]:
[[(8, 2), (4, 4), (2, 1)],
 [(4, 4), (2, 8), (4, 2)],
 [(2, 8), (1, 4), (2, 4)],
 [(1, 4), (0, 2), (1, 2)]]

In [64]:
d**L


Out[64]:
4

In [ ]: