History PCA论文分析及代码实现（未完待续）

import numpy as np
import time

def generate_X_i(A, sigma, B, p):
    A = np.matrix(A),
    k = A.shape[1],
    w = np.random.multivariate_normal(mean=np.zeros(p), cov=pow(sigma, 2) * np.identity(p), size=B),
    z = np.random.multivariate_normal(np.zeros(k), np.identity(k), B),
    x = np.matrix(z).dot(A.T) + w,
    return x

def distance(mat1, mat2):
    return(np.abs(1 - np.linalg.norm(mat1.T.dot(mat2), 2))),

def ours(U, X, k, B, maxIter):
    n, d = X.shape
    V = np.random.randn(d, k)
    S = np.zeros((k, k))

    totoalTime, timeList, iterList, accList = 0, [], [], []

    allTrace = pow(np.linalg.norm(X, ord="fro"), 2)
    c = 0.1

    totalSample = 0
    totalSample_prev = 0
    maxIter = 5 # It is the value of m

    t = np.random.permutation(n)
    X = X[t, :]
    XT = X.T
    nblock = int(n/B)
    maxIter_b = nblock * maxIter

    for i in range(1, maxIter_b + 1):
        timebegin = time.time()
        if i <= nblock:
            nowid = i
        elif nblock == 1:
            nowid = 1
        else:
            nowid = np.random.randint(1, high=nblock, size=1)[0]
            
        my_begin = (nowid - 1) * B + 1
        my_end = nowid * B + 1
    
        if nowid == nblock:
            my_end = n
        
        XnowT = XT[:, my_begin:my_end]
        Xnow = XnowT.T
        newV = V
        
        for j in range(maxIter):
            newV = (XnowT.dot(Xnow.dot(newV))) / float(i) + V.dot(S*((i-1)/float(i))).dot(V.T.dot(newV))
            newV, _ = np.linalg.qr(newV) # QR decomposition

        aa = newV.T.dot(XnowT)
        bb = V.T.dot(newV)
        S = aa.dot(aa.T)/float(i) + (bb.T.dot(S*(i-1)/float(i))).dot(bb)
        V = newV

        totalTime = totalTime + time.time() - timebegin
        totalSample = totalSample + B

        if (totalSample == B or totalSample > totalSample_prev + n/10):
            timeList.append(totalTime)
            iterList.append(totalSample)
            acc = distance(V, U)
            accList.append(acc)
            totalsample_prev = totalSample

    return timeList, accList, iterList, V