In [17]:

    

import numpy as np

x = np.array([ 4., 4., 0, 0, ], dtype=float)
np.outer(x,x)


    

    
Out[17]:





array([4., 4., 0., 0.])

In [15]:

    

def multi_rankupdate(R,x):
    p = np.size(x)
    x = x.T
    
    k=0
    for k in range(p):
        r = np.sqrt(R[k,k]**2 + x[k]**2)
        c = r/R[k,k]
        s = x[k]/R[k,k]
        R[k,k] = r
        R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c
        x[k+1:p]= c*x[k+1:p] - s*R[k, k+1:p]
        
    for k in range(1,p):
        r = R[k,k]
        c = 1.0
        s = 0.
        x[k+1:p] = x[k+1:p]
    return R


    

In [58]:

    

def cholupdate(R,x):
    p = int(np.size(x))
    for k in range(p):
        r = np.sqrt(R[k,k]**2 + x[k]**2)
        c = r/R[k,k]
        s = x[k]/R[k,k]
        R[k,k] = r
        R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c
        x[k+1:p]= c*x[k+1:p] - s*R[k, k+1:p]
    return R


    

In [2]:

    

from sklearn.datasets import make_spd_matrix
A = make_spd_matrix(6,4)
A


    

    
Out[2]:





array([[ 1.7488399 , -1.92816395, -1.39618642, -0.2769755 , -0.52815529,
         0.24624319],
       [-1.92816395,  3.34435465,  2.07258617,  0.4417173 ,  0.84673143,
        -0.35075244],
       [-1.39618642,  2.07258617,  2.1623261 ,  0.25923918,  0.64428255,
        -0.2329581 ],
       [-0.2769755 ,  0.4417173 ,  0.25923918,  0.6147927 ,  0.15112057,
        -0.00692033],
       [-0.52815529,  0.84673143,  0.64428255,  0.15112057,  0.80141217,
        -0.19682322],
       [ 0.24624319, -0.35075244, -0.2329581 , -0.00692033, -0.19682322,
         0.56240547]])

In [118]:

    

L = np.linalg.cholesky(A)
np.dot(L, L.T.conj()) # verify that L * L.H = A


    

    
Out[118]:





array([[ 2.55704765, -0.53197398,  1.3351233 , -1.46980807],
       [-0.53197398,  0.39179579, -0.58107637,  0.14898069],
       [ 1.3351233 , -0.58107637,  1.12295   , -0.70014778],
       [-1.46980807,  0.14898069, -0.70014778,  1.32457697]])

In [132]:

    

beta = np.array([ 1.0, 0, 0, 0 ])
# beta = np.array([ 1.0, 1.0, 1.0, 1.0 ])
A_b = A + np.outer(beta,beta.T)
A_b


    

    
Out[132]:





array([[ 3.55704765, -0.53197398,  1.3351233 , -1.46980807],
       [-0.53197398,  0.39179579, -0.58107637,  0.14898069],
       [ 1.3351233 , -0.58107637,  1.12295   , -0.70014778],
       [-1.46980807,  0.14898069, -0.70014778,  1.32457697]])

In [145]:

    

#function [L] = cholupdate(L, x)
#    n = length(x);
#    for k = 1:n
#        r = sqrt(L(k, k)^2 + x(k)^2);
#        c = r / L(k, k);
#        s = x(k) / L(k, k);
#        L(k, k) = r;
#        L(k+1:n, k) = (L(k+1:n, k) + s * x(k+1:n)) / c;
#        x(k+1:n) = c * x(k+1:n) - s * L(k+1:n, k);
#    end
#end
def cholupdate(Rin,xin):
    R = Rin.copy()
    x = xin.copy()
    p = int(np.size(x))
    for k in range(p):
        r = np.sqrt(R[k,k]**2 + x[k]**2)
        c = r/R[k,k]
        s = x[k]/R[k,k]
        R[k,k] = r
        R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c
        x[k+1:p] = c*x[k+1:p] - s*R[k, k+1:p]
    return R

def cholupdate_allzerosbutone(Rin,xin):
    R = Rin.copy()
    x = xin.copy()
    p = int(np.size(x))
    
    k = 0
    r = np.sqrt(R[k,k]**2 + x[k]**2)
    c = r/R[k,k]
    s = x[k]/R[k,k]
    R[k,k] = r
    R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c
    x[k+1:p] = - s*R[k, k+1:p]
    
    for k in range(1,p):
        print(x[k])
        r = np.sqrt(R[k,k]**2 + x[k]**2)
        c = r/R[k,k]
        s = x[k]/R[k,k]
        R[k,k] = r
        R[k,k+1:p] = (R[k,k+1:p] + s*x[k+1:p])/c
        x[k+1:p] = c*x[k+1:p] - s*R[k, k+1:p]
    return R

L_b = cholupdate_allzerosbutone(L.T, beta).T
np.linalg.norm(np.dot(L_b, L_b.T.conj()) - A_b) # verify that L * L.H = A
# L, L_b


    

    




0.17639088607246817
-0.23947790293899615
0.24510932781116973






    
Out[145]:





2.160125706819257e-15

In [125]:

    

L_b = np.linalg.cholesky(A_b)
L_b


    

    
Out[125]:





array([[ 1.88601369,  0.        ,  0.        ,  0.        ],
       [ 0.24815621,  1.15334916,  0.        ,  0.        ],
       [ 1.23812638,  0.0968266 ,  0.76198273,  0.        ],
       [-0.24910109,  1.0498093 ,  0.66487253,  0.84756744]])

In [ ]: