PSD.py 文件源码

def J(v_X):
    X = vec_to_X(v_X)
    n = X.shape[0]

    # perform scaling
    _i = tuple(range(n))
    X[_i, _i] *= _sqrt2

    Lam, U = np.linalg.eigh(X)
    idx = np.argsort(Lam)
    Lam = Lam[idx]
    U = U[:, idx]
    L = np.diag(Lam)
    L_max = np.maximum(L, 0.0)

    dU_dX, dL_dX = dU_dL_dX(Lam, U)
    dL_max_dX = dL_dX.copy()
    for i, l in enumerate(Lam):
        if l < 0:
            dL_max_dX[i, :, :] = 0
    t1 = dot(U.dot(L_max), np.rollaxis(dU_dX, 1, 0))
    t2 = np.rollaxis(t1, 1, 0)
    t3 = np.rollaxis(dot(multiply_diag(U, dL_max_dX), U.T, (1, 0)), 3, 1)

    idx = np.nonzero(np.tri(n, dtype=np.bool).T)
    W = t1 + t2 + t3

    # rescale jacobian
    W[:, :, _i, _i] *= _sqrt2
    W[_i, _i, :, :] /= _sqrt2

    return W[idx[0], idx[1]][:, idx[0], idx[1]]