logdet_mc.py 文件源码

def update_X(self, X, mu, k=20):
        U, S, VT = svdp(X, k=k)
        P = np.c_[np.ones((k, 1)), 1-S, 1./2./mu-S]
        sigma_star = np.zeros(k)
        for t in range(k):
            p = P[t, :]
            delta = p[1]**2 - 4 * p[0] * p[2]
            if delta <= 0:
                sigma_star[t] = 0.
            else:
                solution = np.roots(p)
                solution = sorted(solution, key=abs)
                solution = np.array(solution)
                if solution[0] * solution[1] <= 0:
                    sigma_star[t] = solution[1]
                elif solution[1] < 0:
                    sigma_star[t] = 0.
                else:
                    f = np.log(1 + solution[1]) + mu * (solution[1] - s[t])**2
                    if f > mu * s[t]**2:
                        sigma_star[t] = 0.
                    else:
                        sigma_star[t] = solution[1]

        sigma_star = np.diag(sigma_star)
        sigma_star = np.dot(np.dot(U, sigma_star), VT)
        return sigma_star