A Closed Form Solution to Robust Subspace Estimation and Clustering中Lemma1中的python实现

最新推荐文章于 2022-04-29 17:00:34 发布

开飞机的小毛驴儿

最新推荐文章于 2022-04-29 17:00:34 发布

阅读量380

点赞数

分类专栏：论文总结

本文链接：https://blog.csdn.net/jzwong/article/details/104079265

版权

论文总结专栏收录该内容

19 篇文章 0 订阅

订阅专栏

用python实现了Lemma1，应该是正确的：

参考代码如下，自己用随机生成的数据做测试，公式得到的最优解和根据C得到的最优解是相同的。

# An implementation of Lemma 1 in the paper
#"A Closed Form Solution to Robust Subspace Estimation and Clustering"
# min||C||*+tau/2||A-AC||

import numpy as np
import pandas as pd
from scipy import linalg
from numpy import *

def svd_with_sigma_as_matrix(input_matrix):
    """
    input_matrix is in 'mat' format
    [u,sigma,v_t]=np.linalg.svd(A), singular values 'sigma' are in a vector format, which means that we can
    not directly use u*sigma*v_t to approximate A. So, we must convert 'sigma' into a full matrix (mxn)
    :param input_matrix [m,n]: matrix that need to decomposed
    :return: u: mxm; sigma:mxn v_t:nxn
    """
    nr, nc = input_matrix.shape
    u, sigma_vector, v_t = np.linalg.svd(input_matrix)
    sigma_full_matrix = np.zeros([nr, nc])
    diag_sigma=np.diag(sigma_vector)
    if nc <= nr:
        sigma_full_matrix[0:nc, :] = diag_sigma[0:nc, :]
    else:
        sigma_full_matrix[:, 0:nr] = diag_sigma[:, 0:nr]

    return u,sigma_full_matrix,sigma_vector,v_t

def closed_form_solution(A,tau):
    U,Sigma_full_matrix,Sigma_vector,V_t=svd_with_sigma_as_matrix(A)
    V=V_t.T

    sqrt_tau=sqrt(tau)

    svp=len(np.flatnonzero(Sigma_vector>1.0/sqrt_tau))

    V1=V[:,0:svp]

    Sigma1=Sigma_vector[0:svp]
    Sigma1_square = multiply(Sigma1, Sigma1)

    Identity_mat=mat(identity(svp))
    print(2*tau*Sigma1_square)
    Phi_I1 = sum(1 - 1.0 / (2 * tau * Sigma1_square))
    Phi_I2=0

    if len(Sigma_vector)-svp>0:
        Sigma2=Sigma_vector[svp:]
        print(len(Sigma2))
        Sigma2_square=diag(multiply(Sigma2,Sigma2))
        Phi_I2+= sum(tau * Sigma2_square / 2.0)

    return V1*(Identity_mat-1.0/tau*diag(1.0/Sigma1_square))*V1.T, Phi_I1+Phi_I2


def main():
    tau=1
    A=mat(random.randint(10,20,size=(500,20)))

    C,optimal_value=closed_form_solution(A,tau)

    print("Low-Rank Matrix:",C)

    print("Rank of C:",np.linalg.matrix_rank(C))

    print("optimal output：",optimal_value)

    [u_c,sigma_c,v_c_t]=linalg.svd(C)
    C_nuclear_norm=sum(sigma_c)

    C_recovery_loss=sum(sum(multiply(A-A*C,A-A*C)))
    print("optimal output manual：", C_nuclear_norm+tau*C_recovery_loss/2.0)



if __name__ == '__main__':
    main()