solver.py 文件源码

def _solve_hessian(G, Y, thY, precon, lambda_min):
    N, T = Y.shape
    # Compute the derivative of the score
    psidY = ne.evaluate('(- thY ** 2 + 1.) / 2.')  # noqa
    # Build the diagonal of the Hessian, a.
    Y_squared = Y ** 2
    if precon == 2:
        a = np.inner(psidY, Y_squared) / float(T)
    elif precon == 1:
        sigma2 = np.mean(Y_squared, axis=1)
        psidY_mean = np.mean(psidY, axis=1)
        a = psidY_mean[:, None] * sigma2[None, :]
        diagonal_term = np.mean(Y_squared * psidY) + 1.
        a[np.diag_indices_from(a)] = diagonal_term
    else:
        raise ValueError('precon should be 1 or 2')
    # Compute the eigenvalues of the Hessian
    eigenvalues = 0.5 * (a + a.T - np.sqrt((a - a.T) ** 2 + 4.))
    # Regularize
    problematic_locs = eigenvalues < lambda_min
    np.fill_diagonal(problematic_locs, False)
    i_pb, j_pb = np.where(problematic_locs)
    a[i_pb, j_pb] += lambda_min - eigenvalues[i_pb, j_pb]
    # Invert the transform
    return (G * a.T - G.T) / (a * a.T - 1.)