GP_multi_classification.py 文件源码

def model_training(K, y, C, n):
    """
    train the model
    :param K: kernel matrix
    :param y: labels which has the same length as all functions f
    :param C: num of categories
    :param n: num of training data
    :return:
    """
    tolerance = 0.0001
    step_size = 0.00001
    s = 0.0005
    # initialization
    f = np.zeros((C*n,)) # initialize f=0(unbiased) which is an constant=0 function and means no GP prior in this case
    block_K = [K[i * n:(i + 1) * n, i * n:(i + 1) * n] for i in range(K.shape[0] / n)]

    # Newton iteration
    for j in range(100):
        pi_vector, pi_matrix = compute_pi(f, C, n)
        D = np.zeros((C * n, C * n))
        np.fill_diagonal(D, pi_vector)
        savetxt_compact('D.txt',D)
        block_D = [D[i * n:(i + 1) * n, i * n:(i + 1) * n] for i in range(D.shape[0] / n)]
        savetxt_compact('block_D0.txt', block_D[0])
        savetxt_compact('block_D2.txt', block_D[2])
        E_c_sum = np.zeros((n, n))
        for c in range(C):
            L = np.linalg.cholesky(np.eye(n) + np.dot(np.sqrt(block_D[c]), np.dot(block_K[c], np.sqrt(block_D[c]))))
            L_inv = np.linalg.inv(L)
            E_c_part = np.dot(np.sqrt(block_D[c]), np.dot(L_inv.T, np.dot(L_inv, np.sqrt(block_D[c]))))
            # create a block diagonal matrix E
            if c == 0:
                E = E_c_part
            else:
                E = block_diag(E, E_c_part)
            E_c_sum += E_c_part
        L_whole = np.linalg.cholesky(np.eye(C*n) + np.dot(np.sqrt(D), np.dot(K, np.sqrt(D))))
        L_whole_inv = np.linalg.inv(L_whole)
        E = np.dot(np.sqrt(D), np.dot(L_whole_inv.T, np.dot(L_whole_inv, np.sqrt(D))))
        #E = np.dot(np.sqrt(D), np.dot(np.linalg.inv(np.eye(C*n) + np.dot(np.sqrt(D), np.dot(K, np.sqrt(D)))), np.sqrt(D)))
        R = np.dot(np.linalg.inv(D), pi_matrix)
        #M = np.linalg.cholesky(E_c_sum)
        M = np.linalg.cholesky(np.dot(R.T, np.dot(E, R)))
        W = D - np.dot(pi_matrix, pi_matrix.T)
        L_K = np.linalg.cholesky(s * np.eye(C*n) + K)
        L_K_inv = np.linalg.inv(L_K)

        b = np.dot((1-step_size) * np.dot(L_K_inv.T, L_K_inv) + W, f) + y - pi_vector
        c = np.dot(E, np.dot(K, b))
        M_inv = np.linalg.inv(M)
        a = b - c + np.dot(E, np.dot(R, np.dot(M_inv.T, np.dot(M_inv, np.dot(R.T, c)))))
        f_new = np.dot(K, a)

        error = np.sqrt(np.sum((f_new - f) ** 2))
        f = f_new
        print `j + 1` + "th iteration, error:" + `error`
        if error <= tolerance:
            print "The function has already converged after " + `j + 1` + " iterations!"
            print "The error is " + `error`
            print "training end!"
            break