diff --git a/README.md b/README.md index 99cc6f8..055732a 100644 --- a/README.md +++ b/README.md @@ -72,12 +72,16 @@ List of algorithms - **LC-KSVD** (Label Consistent K-SVD) - Z. Jiang, Z. Lin, L. S. Davis, "[Learning a discriminative dictionary for sparse coding via label consistent K-SVD](http://ieeexplore.ieee.org/abstract/document/5995354/)," IEEE Conference on Computer Vision and Pattern Recognition (CVPR2011), 2011. - Z. Jiang, Z. Lin, L. S. Davis, "[Label consistent K-SVD: learning A discriminative dictionary for recognition](http://ieeexplore.ieee.org/document/6516503/)," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.35, no.11, pp.2651-2664, 2013. + - **FDDL** (Fisher Discriminative Dictionary Learning) + - M. Yang, L. Zhang, X. Feng, and D. Zhang, "[Fisher discrimination dictionary learning for sparse representation](http://ieeexplore.ieee.org/document/6126286/)," IEEE International Conference on Computer Vision (ICCV), 2011. - **Geometry-aware** - **R-KSRC (Stein kernel)** (Riemannian kernelized sparse representation classification) - M. Harandi, R. Hartley, B. Lovell and C. Sanderson, "[Sparse coding on symmetric positive definite manifolds using bregman divergences](http://ieeexplore.ieee.org/document/7024121/)," IEEE Transactions on Neural Networks and Learning Systems (TNNLS), 2016. - M. Harandi, C. Sanderson, R. Hartley and B. Lovell, "[Sparse coding and dictionary learning for symmetric positive definite matrices: a kernel approach](https://drive.google.com/uc?export=download&id=0B9_PW9TCpxT0eW00U1FVd0xaSmM)," European Conference on Computer Vision (ECCV), 2012. -**R-KSRC (Log-Euclidean kernel)** (Riemannian kernelized sparse representation classification) - P. Li, Q. Wang, W. Zuo, and L. Zhang, "[Log-Euclidean kernels for sparse representation and dictionary learning](http://ieeexplore.ieee.org/document/6751309/)," IEEE International Conference on Computer Vision (ICCV), 2013. + - S. Jayasumana, R. Hartley, M. Salzmann, H. Li, and M. Harandi, "[Kernel methods on the Riemannian manifold of symmetric positive definite matrices](http://ieeexplore.ieee.org/document/6618861/)," IEEE Conference on Computer Vision and Pattern Recognition (CVPR2013), 2013. + - S. Jayasumana, R. Hartley, M. Salzmann, H. Li, and M. Harandi, "[Kernel methods on the Riemannian manifold with Gaussian RBF Kernels](http://ieeexplore.ieee.org/document/7063231/)," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol.37, no.12, 2015. - [Reference] **R-KSRC (Deta-dependent kernel)** [Not included in this package] - Y. Wu, Y. Jia, P. Li, J. Zhang, and J. Yuan, "[Manifold kernel sparse representation of symmetric positive definite matrices and its applications](http://ieeexplore.ieee.org/document/7145428/)," IEEE Transactions on Image Processing, vol.24, no.11, 2015. @@ -199,10 +203,11 @@ Third party tools, libraries, and packages. - [KSVDBox](http://www.cs.technion.ac.il/~ronrubin/Software/ksvdbox13.zip) is used for K-SVD algorithm. - [SPAMS](http://spams-devel.gforge.inria.fr/downloads.html) is used for various lasso problems. - [LC-KSVD](https://www.umiacs.umd.edu/~zhuolin/projectlcksvd.html). + - [FDDL](http://www4.comp.polyu.edu.hk/~cslzhang/code/FDDL.zip). - [RSR](https://drive.google.com/uc?export=download&id=0B9_PW9TCpxT0ZVpGRDNLX3NCbXc). - [Learning Discriminative Stein Kernel for SPD Matrices and Its Applications](https://github.com/seuzjj/DSK/archive/master.zip). - [SDR-SLR](http://www3.ntu.edu.sg/home/EXDJiang/CodesPAMI2015.zip). - - [*R-KSRC (Log-Euclidean kernel)](http://www4.comp.polyu.edu.hk/~cslzhang/LogEKernel_Project/ICCV_LogEKernel_Code.zip). + - [R-KSRC (Log-Euclidean kernel)](http://www4.comp.polyu.edu.hk/~cslzhang/LogEKernel_Project/ICCV_LogEKernel_Code.zip). - DERLR. - [JACOBI_EIGENVALUE](https://people.sc.fsu.edu/~jburkardt/m_src/jacobi_eigenvalue/jacobi_eigenvalue.html) is a MATLAB library which computes the eigenvalues and eigenvectors of a real symmetric matrix. - [NMFLibrary](https://github.com/hiroyuki-kasai/NMFLibrary) is for [NMF](https://en.wikipedia.org/wiki/Non-negative_matrix_factorization). diff --git a/lib/FDDL/Readme.txt b/lib/FDDL/Readme.txt new file mode 100755 index 0000000..e2cbb30 --- /dev/null +++ b/lib/FDDL/Readme.txt @@ -0,0 +1,39 @@ +% ======================================================================== +% Fisher Discriminative Dictionary Learning (FDDL), Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% The code is for the paper: + +% M. Yang, L. Zhang, X. Feng and D. Zhang, +% ˇ§Fisher Discrimination Dictionary Learning for Sparse Representation,ˇ¨ in ICCV 2011. + +% ---------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- + +demo.m Face recognition demo on AR database with 300-d Eigenface feature + +utilier : folder of FDDL functions, including +Eigenface_f: function of computing PCA Projection Matrix +FDDL: main function of FDDL +FDDL_Class_Energy: function of computing energy of certain class +FDDL_FDL_Energy: function of computing energy of all classes +FDDL_Gradient_Comp: function of computing coding model's gradient +FDDL_INIC: function of initializing representation coef +FDDL_INID: function of initializing dictionary +FDDL_SpaCoef: function of computing coding coefficient +FDDL_UpdateDi: function of updating dictioary +IPM_SC: sparse coding function +soft: soft threholding function + +%------------------------------------------------------------------------- +Contact: {csmyang,cslzhang}@comp.polyu.edu.hk \ No newline at end of file diff --git a/lib/FDDL/demo.m b/lib/FDDL/demo.m new file mode 100755 index 0000000..cfa0588 --- /dev/null +++ b/lib/FDDL/demo.m @@ -0,0 +1,57 @@ +close all; +clear all; +clc; + +addpath([cd '/utilies']); +load(['AR_EigenFace']); + +%%%%%%%%%%%%%%%%%%%%%%%% +%FDDL parameter +%%%%%%%%%%%%%%%%%%%%%%%% +opts.nClass = 100; +opts.wayInit = 'PCA'; +opts.lambda1 = 0.005; +opts.lambda2 = 0.05; +opts.nIter = 15; +opts.show = true; +[Dict,Drls,CoefM,CMlabel] = FDDL(tr_dat,trls,opts); + +%%%%%%%%%%%%%%%%%%%%%%%% +% Sparse Classification +%%%%%%%%%%%%%%%%%%%%%%%% +lambda = 0.005; +nClass = opts.nClass; +weight = 0.5; + +td1_ipts.D = Dict; +td1_ipts.tau1 = lambda; +if size(td1_ipts.D,1)>=size(td1_ipts.D,2) + td1_par.eigenv = eigs(td1_ipts.D'*td1_ipts.D,1); +else + td1_par.eigenv = eigs(td1_ipts.D*td1_ipts.D',1); +end + +ID = []; +for indTest = 1:size(tt_dat,2) + fprintf(['Totalnum:' num2str(size(tt_dat,2)) 'Nowprocess:' num2str(indTest) '\n']); + td1_ipts.y = tt_dat(:,indTest); + [opts] = IPM_SC(td1_ipts,td1_par); + s = opts.x; + + for indClass = 1:nClass + temp_s = zeros(size(s)); + temp_s(indClass==Drls) = s(indClass==Drls); + zz = tt_dat(:,indTest)-td1_ipts.D*temp_s; + gap(indClass) = zz(:)'*zz(:); + + mean_coef_c = CoefM(:,indClass); + gCoef3(indClass) = norm(s-mean_coef_c,2)^2; + end + + wgap3 = gap + weight*gCoef3; + index3 = find(wgap3==min(wgap3)); + id3 = index3(1); + ID = [ID id3]; +end + +fprintf('%s%8f\n','reco_rate = ',sum(ID==ttls)/(length(ttls))); \ No newline at end of file diff --git a/lib/FDDL/utilies/Eigenface_f.m b/lib/FDDL/utilies/Eigenface_f.m new file mode 100755 index 0000000..9304cb6 --- /dev/null +++ b/lib/FDDL/utilies/Eigenface_f.m @@ -0,0 +1,54 @@ +function [disc_set,disc_value,Mean_Image]=Eigenface_f(Train_SET,Eigen_NUM) + +% the magnitude of eigenvalues of this function is corrected right !!!!!!!!! +% Centralized PCA +[NN,Train_NUM]=size(Train_SET); + +if NN<=Train_NUM % for small sample size case + + Mean_Image=mean(Train_SET,2); + Train_SET=Train_SET-Mean_Image*ones(1,Train_NUM); + R=Train_SET*Train_SET'/(Train_NUM-1); + + [V,S]=Find_K_Max_Eigen(R,Eigen_NUM); + disc_value=S; + disc_set=V; + +else % for small sample size case + + Mean_Image=mean(Train_SET,2); + Train_SET=Train_SET-Mean_Image*ones(1,Train_NUM); + + R=Train_SET'*Train_SET/(Train_NUM-1); + + [V,S]=Find_K_Max_Eigen(R,Eigen_NUM); + disc_value=S; + disc_set=zeros(NN,Eigen_NUM); + + Train_SET=Train_SET/sqrt(Train_NUM-1); + for k=1:Eigen_NUM + disc_set(:,k)=(1/sqrt(disc_value(k)))*Train_SET*V(:,k); + end + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [Eigen_Vector,Eigen_Value]=Find_K_Max_Eigen(Matrix,Eigen_NUM) + +[NN,NN]=size(Matrix); +[V,S]=eig(Matrix); %Note this is equivalent to; [V,S]=eig(St,SL); also equivalent to [V,S]=eig(Sn,St); % + +S=diag(S); +[S,index]=sort(S); + +Eigen_Vector=zeros(NN,Eigen_NUM); +Eigen_Value=zeros(1,Eigen_NUM); + +p=NN; +for t=1:Eigen_NUM + Eigen_Vector(:,t)=V(:,index(p)); + Eigen_Value(t)=S(p); + p=p-1; +end diff --git a/lib/FDDL/utilies/FDDL.m b/lib/FDDL/utilies/FDDL.m new file mode 100755 index 0000000..a8dcd30 --- /dev/null +++ b/lib/FDDL/utilies/FDDL.m @@ -0,0 +1,170 @@ +function [Dict,Drls,CoefM,CMlabel] = FDDL(TrainDat,TrainLabel,opts) +% ======================================================================== +% Fisher Discriminative Dictionary Learning (FDDL), Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ---------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- +% +% This is an implementation of the algorithm for learning the +% Fisher Discriminative Dictionary from a labeled training data +% +% Please refer to the following paper +% +% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination +% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on +% Computer Vision, 2011. +% +%---------------------------------------------------------------------- +% +%Input : (1) TrainDat: the training data matrix. +% Each column is a training sample +% (2) TrainDabel: the training data labels +% (3) opts : the struture of parameters +% .nClass the number of classes +% .wayInit the way to initialize the dictionary +% .lambda1 the parameter of l1-norm energy of coefficient +% .lambda2 the parameter of l2-norm of Fisher Discriminative +% coefficient term +% .nIter the number of FDDL's iteration +% .show sign value of showing the gap sequence +% +%Output: (1) Dict: the learnt dictionary via FDDL +% (2) Drls: the labels of learnt dictionary's columns +% (2) CoefM: Mean Coefficient Matrix. Each column is a mean coef +% vector +% (3) CMlabel: the labels of CoefM's columns. +% +%----------------------------------------------------------------------- +% +%Usage: +%Given a training data, including TrainDat and TrainLabel, and the +%parameters, opts. +% +%[Dict,CoefM,CMlabel] = FDDL(TrainDat,TrainLabel,opts) +%----------------------------------------------------------------------- +%%%%%%%%%%%%%%%%%% +% normalize energy +%%%%%%%%%%%%%%%%%% +TrainDat = TrainDat*diag(1./sqrt(sum(TrainDat.*TrainDat))); + +%%%%%%%%%%%%%%%%%% +%initialize dict +%%%%%%%%%%%%%%%%%% +Dict_ini = []; +Dlabel_ini = []; +for ci = 1:opts.nClass + cdat = TrainDat(:,TrainLabel==ci); + dict = FDDL_INID(cdat,size(cdat,2),opts.wayInit); + Dict_ini = [Dict_ini dict]; + Dlabel_ini = [Dlabel_ini repmat(ci,[1 size(dict,2)])]; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%initialize coef without between-class scatter +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +ini_par.tau = opts.lambda1; +ini_par.lambda = opts.lambda2; +ini_ipts.D = Dict_ini; +coef = zeros(size(Dict_ini,2),size(TrainDat,2)); +if size(Dict_ini,1)>size(Dict_ini,2) + ini_par.c = 1.05*eigs(Dict_ini'*Dict_ini,1); +else + ini_par.c = 1.05*eigs(Dict_ini*Dict_ini',1); +end +for ci = 1:opts.nClass + fprintf(['Initializing Coef: Class ' num2str(ci) '\n']); + ini_ipts.X = TrainDat(:,TrainLabel==ci); + [ini_opts] = FDDL_INIC (ini_ipts,ini_par); + coef(:,TrainLabel ==ci) = ini_opts.A; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%Main loop of Fisher Discriminative Dictionary Learning +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + Fish_par.dls = Dlabel_ini; + Fish_ipts.D = Dict_ini; + Fish_ipts.trls = TrainLabel; + Fish_par.tau = opts.lambda1; + Fish_par.lambda2 = opts.lambda2; + + Fish_nit = 1; + drls = Dlabel_ini; + while Fish_nit<=opts.nIter + if size(Fish_ipts.D,1)>size(Fish_ipts.D,2) + Fish_par.c = 1.05*eigs(Fish_ipts.D'*Fish_ipts.D,1); + else + Fish_par.c = 1.05*eigs(Fish_ipts.D*Fish_ipts.D',1); + end + %------------------------- + %updating the coefficient + %------------------------- + for ci = 1:opts.nClass + fprintf(['Updating coefficients, class: ' num2str(ci) '\n']) + Fish_ipts.X = TrainDat(:,TrainLabel==ci); + Fish_ipts.A = coef; + Fish_par.index = ci; + [Copts] = FDDL_SpaCoef (Fish_ipts,Fish_par); + coef(:,TrainLabel==ci) = Copts.A; + CMlabel(ci) = ci; + CoefM(:,ci) = mean(Copts.A,2); + end + [GAP_coding(Fish_nit)] = FDDL_FDL_Energy(TrainDat,coef,opts.nClass,Fish_par,Fish_ipts); + + %------------------------ + %updating the dictionary + %------------------------ + for ci = 1:opts.nClass + fprintf(['Updating dictionary, class: ' num2str(ci) '\n']); + [Fish_ipts.D(:,drls==ci),Delt(ci).delet]= FDDL_UpdateDi (TrainDat,coef,... + ci,opts.nClass,Fish_ipts,Fish_par); + end + [GAP_dict(Fish_nit)] = FDDL_FDL_Energy(TrainDat,coef,opts.nClass,Fish_par,Fish_ipts); + + newD = []; newdrls = []; newcoef = []; + for ci = 1:opts.nClass + delet = Delt(ci).delet; + if isempty(delet) + classD = Fish_ipts.D(:,drls==ci); + newD = [newD classD]; + newdrls = [newdrls repmat(ci,[1 size(classD,2)])]; + newcoef = [newcoef; coef(drls==ci,:)]; + else + temp = Fish_ipts.D(:,drls==ci); + temp_coef = coef(drls==ci,:); + for temp_i = 1:size(temp,2) + if sum(delet==temp_i)==0 + newD = [newD temp(:,temp_i)]; + newdrls = [newdrls ci]; + newcoef = [newcoef;temp_coef(temp_i,:)]; + end + end + end + end + + Fish_ipts.D = newD; + coef = newcoef; + drls = newdrls; + Fish_par.dls = drls; + + Fish_nit = Fish_nit +1; + end + + Dict = Fish_ipts.D; + Drls = drls; + + if opts.show + subplot(1,2,1);plot(GAP_coding,'-*');title('GAP_coding'); + subplot(1,2,2);plot(GAP_dict,'-o');title('GAP_dict'); + end +return; \ No newline at end of file diff --git a/lib/FDDL/utilies/FDDL_Class_Energy.m b/lib/FDDL/utilies/FDDL_Class_Energy.m new file mode 100755 index 0000000..a9f520a --- /dev/null +++ b/lib/FDDL/utilies/FDDL_Class_Energy.m @@ -0,0 +1,60 @@ +function [gap] = FDDL_Class_Energy(Ai,D,Xi,Xa,drls,trls,index,lambda1,lambda2,lambda3,lambda4,classn,tau2,tau3) +% ======================================================================== +% Class energy computation of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% +% Input : (1) Ai : the data matrix of this class +% (2) D : the whole dictionary +% (3) Xi: the coefficient matrix of this class +% (4) Xa: the coefficient matrix of the whole class +% (5) drls: labels of dictionary's column +% (6) trls: labels of training samples +% (7) index: label of class being processed +% (8) lambda1 : parameter of l1-norm energy of coefficient +% (9) lambda2 : parameter of within-class scatter +% (10) lambda3 : parameter of between-class scatter +% (11) lambda4: parameter of l2-norm energy of coefficient +% (12) classn: the number of class +% (13) tau2: parameter of ||A_i-D_iX_i^i||_F^2 in fidelity term +% (14) tau3: parameter of ||D_jX_i^j||_F^2 in fidelity term +% +% Outputs : (1) gap : the total energy of some class +% +%------------------------------------------------------------------------ + +gap3 = 0; +gap4 = 0; +GAP1 = norm((Ai-D*Xi),'fro')^2;% only for one class, no effect +GAP2 = lambda1*sum(abs(Xi(:)));% + +Xa(:,trls==index) = Xi; +tem_XA = Xa; +mean_xa = mean(tem_XA,2); + + n_c = size(Xi,2); + for i_c = 1:classn + t_X_ic = tem_XA(:,trls==i_c); + n_ic = size(t_X_ic,2); + mean_xic = mean(t_X_ic,2); +% gap3 = gap3+norm(t_X_ic-repmat(mean(t_X_ic,2),[1 size(t_X_ic,2)]),'fro')^2; + gap4 = gap4+n_ic*(mean_xic-mean_xa)'*(mean_xic-mean_xa); + end + + gap3 = norm(Xi-repmat(mean(Xi,2),[1 n_c]),'fro')^2; + + GAP3 = lambda2*gap3-lambda3*gap4; + + GAP4 = lambda4*norm(Xi,'fro')^2;% only for one class, no effect + + GAP5 = (tau2)*norm((Ai-D(:,drls==index)*Xi(drls==index,:)),'fro')^2;% only for one class, no effect + GAP6 = 0; + for i = 1:classn + if i~=index + GAP6 = GAP6+tau3*norm(D(:,drls==i)*Xi(drls==i,:),'fro')^2; + end + end + + gap = GAP1+GAP2+GAP3+GAP5+GAP6+GAP4; \ No newline at end of file diff --git a/lib/FDDL/utilies/FDDL_FDL_Energy.m b/lib/FDDL/utilies/FDDL_FDL_Energy.m new file mode 100755 index 0000000..2de0ee3 --- /dev/null +++ b/lib/FDDL/utilies/FDDL_FDL_Energy.m @@ -0,0 +1,62 @@ +function [gap] = FDDL_FDL_Energy(Aa,Xa,classn,Fish_par,Fish_ipts) +% ======================================================================== +% Total energy computation of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% +% Input : (1) Aa : the data matrix of all class +% (2) Xa: the coefficient matrix of the whole class +% (3) classn: the number of class +% (4) Fish_par +% .dls labels of dictionary's column +% .tau parameter of l1-norm energy of coefficient +% .lambda2 parameter of within-class scatter +% (5) Fish_ipts +% .D The dictioanry +% .trls labels of training samples +% +% Outputs : (1) gap : the total energy +% +%------------------------------------------------------------------------ + D = Fish_ipts.D; + drls = Fish_par.dls; + trls = Fish_ipts.trls; + lambda1 = Fish_par.tau; + lambda2 = Fish_par.lambda2; + lambda3 = lambda2; + lambda4 = lambda2; + tau2 = 1; + tau3 = 1; + + + gap3 = 0; + gap4 = 0; + GAP1 = norm((Aa-D*Xa),'fro')^2;% + GAP2 = lambda1*sum(abs(Xa(:)));% + tem_XA = Xa; + for i_c = 1:classn + t_X_ic = tem_XA(:,trls==i_c); + gap3 = gap3+norm(t_X_ic-repmat(mean(t_X_ic,2),[1 size(t_X_ic,2)]),'fro')^2; + gap4 = gap4+size(t_X_ic,2)*(mean(t_X_ic,2)-mean(tem_XA,2))'*(mean(t_X_ic,2)-mean(tem_XA,2)); + end + + GAP3 = lambda2*gap3-lambda3*gap4; + + GAP4 = lambda4*norm(Xa,'fro')^2;% + + GAP6 = 0; + GAP5 = 0; + for i = 1:classn + Ai = Aa(:,trls==i); + Xi = Xa(:,trls==i); + GAP5 = GAP5 + (tau2)*norm(Ai-D(:,drls==i)*Xi(drls==i,:),'fro')^2;% only for one class, no effect + for j = 1:classn + if j~=i + GAP6 = GAP6+tau3*norm(D(:,drls==j)*Xi(drls==j,:),'fro')^2; + end + end + end + + gap = GAP1+GAP2+GAP3+GAP5+GAP6+GAP4; \ No newline at end of file diff --git a/lib/FDDL/utilies/FDDL_Gradient_Comp.m b/lib/FDDL/utilies/FDDL_Gradient_Comp.m new file mode 100755 index 0000000..991dd33 --- /dev/null +++ b/lib/FDDL/utilies/FDDL_Gradient_Comp.m @@ -0,0 +1,87 @@ +function grad = FDDL_Gradient_Comp(Xi,Xa,classn,index,lambda2,lambda3,lambda4,tau2,... + tau3,trls,drls,newpar,BAI,CJ) +% ======================================================================== +% IPM's Gradient computation of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% +% Input : (1) Xi: the coefficient matrix of this class +% (2) Xa: the coefficient matrix of the whole class +% (3) classn: the number of class +% (4) index: label of class being processed +% (5) lambda2 : parameter of within-class scatter +% (6) lambda3 : parameter of between-class scatter +% (7) lambda4: parameter of l2-norm energy of coefficient +% (8) tau2: parameter of ||A_i-D_iX_i^i||_F^2 in fidelity term +% (9) tau3: parameter of ||D_jX_i^j||_F^2 in fidelity term +% (10) trls: labels of training samples +% (11) drls: labels of dictionary's column +% (12) newpar, BAI, and CJ: the precomputed data +% +% Outputs : (1) grad : the gradient vector of coding model +% +%------------------------------------------------------------------------ + +n_d = newpar.n_d; % the sample number of i-th training data +B_line_i = newpar.B_line_i; +C_j = newpar.C_j; +C_line = newpar.C_line; +DD = newpar.DD; +DAi = newpar.DAi; +Di0Di0 = newpar.Di0Di0; +Di0Ai = newpar.Di0Ai; +BiBi = newpar.BiBi; +BaiBai = newpar.BaiBai; +BaiGxi = newpar.BaiGxi; +CjCj = newpar.CjCj; +m = newpar.m; +DoiDoi = newpar.DoiDoi; +B_angle_i = newpar.Bai; + +for k = 1:classn +% Z(k).Matrix = Xa(:,trls==k)*B_line_i-Xa*C_line+Xi*C_j; +% Z(k).Matrix = Xa(:,trls==k)*B_angle_i-Xa*C_line+Xi*C_j+Xa(:,trls==k)*C_j; + Z(k).Matrix = Xa(:,trls==k)*BAI(k).M-Xa*C_line+Xi*C_j+Xa(:,trls==k)*CJ(k).M; +end + +XiT = Xi'; + +tem = 2*DD*Xi-2*DAi; +grad1 = tem(:); + +tem = tau3*2*DoiDoi*Xi+(tau2)*2*Di0Di0*Xi-(tau2)*2*Di0Ai; +grad7 = tem(:); + +% grad8 = zeros(size(grad7)); +% for i =1:classn +% Di0 = D; +% Di0(:,drls~=i)=0; +% tem = -(tau3*2*Di0'*Di0*Xi-tau3*2*Di0'*Ai); +% grad8 = grad8+tem(:); +% end + +tem = 2*lambda2*BiBi*XiT; +grad4 = tem(:); + +tem = 2*lambda4*XiT; +grad9 = tem(:); + +tem = -lambda3*(2*BaiBai*XiT-2*BaiGxi); +grad5 = tem(:); + +grad6 = zeros(size(grad5)); +for k = 1:classn + temz = Z(k).Matrix; + if k~=index + tem = -lambda3*(2*CjCj*XiT-2*C_j*temz'); + grad6 = grad6+tem(:); + end +end + +grad456 = reshape(grad4+grad5+grad6+grad9,[n_d m])'; +grad456 = grad456(:); + +% grad = grad1+grad456+grad7+grad8; +grad = grad1+grad456+grad7; diff --git a/lib/FDDL/utilies/FDDL_INIC.m b/lib/FDDL/utilies/FDDL_INIC.m new file mode 100755 index 0000000..8aa7385 --- /dev/null +++ b/lib/FDDL/utilies/FDDL_INIC.m @@ -0,0 +1,234 @@ +function [opts] = FDDL_INIC (ipts,par) +% ======================================================================== +% Coefficient Initialization of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- +% +% This is an implementation of the algorithm for initializing the +% Coefficient matrix of FDDL +% +% Please refer to the following paper +% +% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination +% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on +% Computer Vision, 2011. +% L. Rosasco, A. Verri, M. Santoro, S. Mosci, and S. Villa. Iterative +% Projection Methods for Structured Sparsity Regularization. MIT Technical +% Reports, MIT-CSAIL-TR-2009-050,CBCL-282, 2009. +% J. Bioucas-Dias, M. Figueiredo, ?A new TwIST: two-step iterative shrinkage +% /thresholding algorithms for image restoration?, IEEE Transactions on +% Image Processing, December 2007. +%---------------------------------------------------------------------- +% +% Inputs : (1) ipts : the structre of input data +% .D the dictionary +% .X the training data +% .last_coef the coef in the last iteration +% (2) par : the struture of input parameters +% .tau the parameter of sparse constraint of coef +% .lambda the parameter of within-class scatter +% +% Outputs: (1) opts : the structure of output data +% .A the coefficient matrix +% .ert the total energy sequence +% +%--------------------------------------------------------------------- + +par.initM = 'zero'; % initialization method +par.nIter = 200; % maximal iteration number +par.isshow = true; % +par.twist = true; % 'true': use twist +par.citeT = 1e-6; % stop criterion +par.cT = 1e+10; % stop criterion + +m = size(ipts.D,2); +n = size(ipts.X,2); + +switch lower(par.initM) + case {'zero'} + A = zeros(m,n); + case {'transpose'} + A = ipts.D'*ipts.X; + case {'pinv'} + A = pinv(ipts.D)*ipts.X; + case {'last'} + A = ipts.last_coef; + otherwise + error('Nonknown method!'); +end + +D = ipts.D; +X = ipts.X; +tau = par.tau; +lambda = par.lambda; +nIter = par.nIter; +c = par.c; +sigma = c; +tau1 = tau/2; +B = eye(n)-ones(n,n)/n; + +At_pref = A(:); +At_now = A(:); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%TWIST parameter +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for_ever = 1; +IST_iters = 0; +TwIST_iters = 0; +sparse = 1; +verbose = 1; +enforceMonotone = 1; +lam1 = 1e-4; %default minimal eigenvalues +lamN = 1; %default maximal eigenvalues +rho0 = (1-lam1/lamN)/(1+lam1/lamN); +alpha = 2/(1+sqrt(1-rho0^2)); %default,user can set +beta = alpha*2/(lam1+lamN); %default,user can set + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%main loop +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +xm2 = At_pref; +xm1 = At_pref; + +A = reshape(At_pref,[m,n]); +gap1 = norm((X-D*A),'fro')^2; +if n==1 + gap2 = norm(A*B,2)^2; +else + gap2 = norm(A*B,'fro')^2; +end +gap3 = sum(abs(A(:))); +prev_f = gap1+2*tau1*gap3+lambda*gap2; + +for n_it = 2 : nIter; + A = reshape(At_now,[m,n]); + gap1 = norm((X-D*A),'fro')^2; + if n==1 + gap2 = norm(A*B,2)^2; + else + gap2 = norm(A*B,'fro')^2; + end + gap3 = sum(abs(A(:))); + ert(n_it-1) = gap1+2*tau1*gap3+lambda*gap2; + +% fprintf('Iteration:%f Total gap:%f\n',n_it,ert(n_it-1)); + + while for_ever + % IPM estimate + v1 = []; + for i = 1:n + A = reshape(xm1,[m,n]); + tem1 = X(:,i)-D*A(:,i); + tem2 = D'*tem1; + v1 = [v1;tem2]; + end + A = reshape(xm1,[m,n])'; + v2_2 = []; + for i = 1:m + tem1 = B*A(:,i); + v2_2 = [v2_2;tem1]; + end + v2_3 = reshape(v2_2,[n m])'; + v2 = v2_3(:); + + v = xm1+(v1-lambda*v2)/sigma; + x_temp = soft(v,tau1/sigma); + + if (IST_iters >= 2) | ( TwIST_iters ~= 0) + % set to zero the past when the present is zero + % suitable for sparse inducing priors + if sparse + mask = (x_temp ~= 0); + xm1 = xm1.* mask; + xm2 = xm2.* mask; + end + % two-step iteration + xm2 = (alpha-beta)*xm1 + (1-alpha)*xm2 + beta*x_temp; + % compute residual + + A = reshape(xm2,[m,n]); + gap1 = norm((X-D*A),'fro')^2; + if n==1 + gap2 = norm(A*B,2)^2; + else + gap2 = norm(A*B,'fro')^2; + end + gap3 = sum(abs(A(:))); + f = gap1+2*tau1*gap3+lambda*gap2; + + if (f > prev_f) & (enforceMonotone) + TwIST_iters = 0; % do a IST iteration if monotonocity fails + else + TwIST_iters = TwIST_iters+1; % TwIST iterations + IST_iters = 0; + x_temp = xm2; + if mod(TwIST_iters,10000) ==0 + c = 0.9*c; + sigma= c; + end + break; % break loop while + end + else + A = reshape(x_temp,[m,n]); + gap1 = norm((X-D*A),'fro')^2; + if n==1 + gap2 = norm(A*B,2)^2; + else + gap2 = norm(A*B,'fro')^2; + end + gap3 = sum(abs(A(:))); + f = gap1+2*tau1*gap3+lambda*gap2; + + if f > prev_f + % if monotonicity fails here is because + % max eig (A'A) > 1. Thus, we increase our guess + % of max_svs + c = 2*c; + sigma = c; + if verbose +% fprintf('Incrementing c=%2.2e\n',c); + end + if c > par.cT + break; % break loop while + end + IST_iters = 0; + TwIST_iters = 0; + else + TwIST_iters = TwIST_iters + 1; + break; % break loop while + end + end + end + + citerion = abs(f-prev_f)/prev_f; + if citerion < par.citeT | c > par.cT +% fprintf('Stop!\n c=%2.2e\n citerion=%2.2e\n',c,citerion); + break; + end + + xm2 = xm1; + xm1 = x_temp; + At_pref = At_now; + At_now = x_temp; + prev_f = f; + +end + +opts.A = reshape(At_now,[m,n]); +opts.ert = ert; +if par.isshow + plot(ert,'r-'); +end \ No newline at end of file diff --git a/lib/FDDL/utilies/FDDL_INID.m b/lib/FDDL/utilies/FDDL_INID.m new file mode 100755 index 0000000..131c6e7 --- /dev/null +++ b/lib/FDDL/utilies/FDDL_INID.m @@ -0,0 +1,30 @@ +function D = FDDL_INID(data,nCol,wayInit) +% ======================================================================== +% Dictionary Initialization of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% +% Input : (1) data : the data matrix +% (2) nCol : the number of dictioanry's columns +% (3) wayInit: the method to initialize dictionary +% +% Outputs : (1) D : the initialized dictionary +% +%------------------------------------------------------------------------ + +m = size(data,1); + +switch lower(wayInit) + case {'pca'} + [D,disc_value,Mean_Image] = Eigenface_f(data,nCol-1); + D = [D Mean_Image./norm(Mean_Image)]; + case {'random'} + phi = randn(m, nCol); + phinorm = sqrt(sum(phi.*phi, 2)); + D = phi ./ repmat(phinorm, 1, nCol); + otherwise + error{'Unkonw method.'} +end +return; diff --git a/lib/FDDL/utilies/FDDL_SpaCoef.m b/lib/FDDL/utilies/FDDL_SpaCoef.m new file mode 100755 index 0000000..a15208b --- /dev/null +++ b/lib/FDDL/utilies/FDDL_SpaCoef.m @@ -0,0 +1,246 @@ +function [opts] = FDDL_SpaCoef(ipts,par) +% ======================================================================== +% Coefficient updating of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- +% +% This is an implementation of the algorithm for updating the +% Coefficient matrix of FDDL (fix the dictionary) +% +% Please refer to the following paper +% +% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination +% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on +% Computer Vision, 2011. +% L. Rosasco, A. Verri, M. Santoro, S. Mosci, and S. Villa. Iterative +% Projection Methods for Structured Sparsity Regularization. MIT Technical +% Reports, MIT-CSAIL-TR-2009-050,CBCL-282, 2009. +% J. Bioucas-Dias, M. Figueiredo, ?A new TwIST: two-step iterative shrinkage +% /thresholding algorithms for image restoration?, IEEE Transactions on +% Image Processing, December 2007. +%---------------------------------------------------------------------- +% +% Inputs : (1) ipts : the structre of input data +% .D the dictionary +% .X the training data +% .A the coefficient matrix in the last iteration +% .trls the labels of training data +% (2) par : the struture of input parameters +% .tau the parameter of sparse constraint of coef +% .lambda the parameter of within-class scatter +% .dls the labels of dictionary's columns +% .index the label of the class being processed +% +% Outputs: (1) opts : the structure of output data +% .A the coefficient matrix +% .ert the total energy sequence +% +%--------------------------------------------------------------------- + +par.nIter = 200; % maximal iteration number +par.isshow = false; +par.citeT = 1e-6; % stop criterion +par.cT = 1e+10; % stop criterion + +m = size(ipts.D,2); +fish_tau3 = 1; % parameter of ||A_i-D_iX_i^i||_F^2 in fidelity term +fish_tau2 = 1; % parameter of ||D_jX_i^j||_F^2 in fidelity term Eq.(4) +drls = par.dls; +D = ipts.D; +X = ipts.X; +A = ipts.A; +tau = par.tau; +lambda1 = par.tau; +lambda2 = par.lambda2; +lambda3 = par.lambda2; +lambda4 = par.lambda2; +trls = ipts.trls; +classn = length(unique(trls)); +nIter = par.nIter; +c = par.c; +sigma = c; +tau1 = tau/2; +index = par.index; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%TWIST parameter +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for_ever = 1; +IST_iters = 0; +TwIST_iters = 0; +sparse = 1; +verbose = 1; +enforceMonotone = 1; +lam1 = 1e-4; %default minimal eigenvalues +lamN = 1; %default maximal eigenvalues +rho0 = (1-lam1/lamN)/(1+lam1/lamN); +alpha = 2/(1+sqrt(1-rho0^2)); %default,user can set +beta = alpha*2/(lam1+lamN); %default,user can set + + +%%%%%%%%%%%%%%%%%%%%%%%% +%preprocessing +%%%%%%%%%%%%%%%%%%%%%%%% +Ai = X; % the i-th training data +Xa = A; % +Xi = A(:,trls==index); +Xt_now = A(:,trls==index); +% Xi = zeros(size(A(:,trls==index))); +% Xt_now = zeros(size(A(:,trls==index))); + +newpar.n_d = size(Ai,2); % the sample number of i-th training data +n = size(Xa,2); % the total sample number of training data + +for ci = 1:classn + t_n_d = sum(trls==ci); + t_b_line_i = ones(t_n_d,newpar.n_d)./t_n_d; + t_c_j = ones(t_n_d,newpar.n_d)./n; + t_b_angle_i = t_b_line_i-t_c_j; + CJ(ci).M = t_c_j; + BAI(ci).M = t_b_angle_i; +end + +newpar.B_line_i = ones(newpar.n_d,newpar.n_d)./newpar.n_d; +newpar.C_j = ones(newpar.n_d,newpar.n_d)./n; +newpar.CjCj = (newpar.C_j)*(newpar.C_j)'; +newpar.C_line = ones(n,newpar.n_d)./n; +B_i = eye(newpar.n_d,newpar.n_d)-newpar.B_line_i; +newpar.BiBi = B_i*(B_i)'; +B_angle_i = newpar.B_line_i-newpar.C_j; +newpar.Bai = B_angle_i; +newpar.BaiBai = B_angle_i*(B_angle_i)'; +Xo = Xa; +Xo(:,trls==index) = 0; +G_X_i = Xo*newpar.C_line; +newpar.BaiGxi = B_angle_i*(G_X_i)'; +newpar.DD = D'*D; +newpar.DAi = D'*Ai; +Di0 = D; +Di0(:,drls~=index) = 0; +newpar.Di0Di0 = (Di0)'*Di0; +newpar.Di0Ai = (Di0)'*Ai; + +newpar.DoiDoi = zeros(size(D,2)); +for t_i = 1:classn + if t_i ~= index + Doi = D; + Doi(:,drls~=t_i) = 0; + newpar.DoiDoi = newpar.DoiDoi+(Doi)'*Doi; + end +end + +newpar.m = m; % the number of dictionary column atoms + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%main loop +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +Xa(:,trls==index) = Xi; +xm2 = Xi;%A(:,trls==index); +xm1 = Xi;%A(:,trls==index); % now + +[gap] = FDDL_Class_Energy(Ai,D,xm1,Xa,drls,trls,index,... + lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3); +prev_f = gap; +ert(1) = gap; +for n_it = 2 : nIter; + + Xa(:,trls==index) = Xi; + + while for_ever + % IPM estimate + + grad = FDDL_Gradient_Comp(xm1,Xa,classn,index,... + lambda2,lambda3,lambda4,fish_tau2,fish_tau3,trls,drls,newpar,... + BAI,CJ); + + v = xm1(:)-grad./(2*sigma); + tem = soft(v,tau1/sigma); + x_temp = reshape(tem,[size(D,2),size(xm1,2)]); + + if (IST_iters >= 2) | ( TwIST_iters ~= 0) + % set to zero the past when the present is zero + % suitable for sparse inducing priors + if sparse + mask = (x_temp ~= 0); + xm1 = xm1.* mask; + xm2 = xm2.* mask; + end + % two-step iteration + xm2 = (alpha-beta)*xm1 + (1-alpha)*xm2 + beta*x_temp; + % compute residual + [gap] = FDDL_Class_Energy(Ai,D,xm2,Xa,drls,trls,index,... + lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3); + + f = gap; + + if (f > prev_f) & (enforceMonotone) + TwIST_iters = 0; % do a IST iteration if monotonocity fails + else + TwIST_iters = TwIST_iters+1; % TwIST iterations + IST_iters = 0; + x_temp = xm2; + if mod(TwIST_iters,10000) ==0 + c = 0.9*c; + sigma = c; + end + break; % break loop while + end + else + + [gap] = FDDL_Class_Energy(Ai,D,x_temp,Xa,drls,trls,index,... + lambda1,lambda2,lambda3,lambda4,classn,fish_tau2,fish_tau3); + + f = gap; + + if f > prev_f + % if monotonicity fails here is because + % max eig (A'A) > 1. Thus, we increase our guess + % of max_svs + c = 2*c; + sigma = c; + if verbose +% fprintf('Incrementing c=%2.2e\n',c); + end + if c > par.cT + break; % break loop while + end + IST_iters = 0; + TwIST_iters = 0; + else + TwIST_iters = TwIST_iters + 1; + break; % break loop while + end + end + + end + + citerion = abs(f-prev_f)/abs(prev_f); + if citerion < par.citeT | c > par.cT +% fprintf('Stop!\n c=%2.2e\n citerion=%2.2e\n',c,citerion); + break; + end + + xm2 = xm1; + xm1 = x_temp; + Xt_now = x_temp; + Xi = Xt_now; + prev_f = f; + ert(n_it) = f; +% fprintf('Iteration:%f Total gap:%f\n',n_it,ert(n_it-1)); + end + + +opts.A = Xt_now; +opts.ert = ert; \ No newline at end of file diff --git a/lib/FDDL/utilies/FDDL_UpdateDi.m b/lib/FDDL/utilies/FDDL_UpdateDi.m new file mode 100755 index 0000000..71e6ef3 --- /dev/null +++ b/lib/FDDL/utilies/FDDL_UpdateDi.m @@ -0,0 +1,95 @@ +function [Di_new,delet] = FDDL_UpdateDi (A,X,index,classn,Fish_ipts,Fish_par) +% ======================================================================== +% Dictionary updating of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- +% +% This is an implementation of the algorithm for updating the +% dictionary of FDDL (fix the coefficient) +% +% Please refer to the following paper +% +% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination +% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on +% Computer Vision, 2011. +% Meng Yang, Lei Zhang, Jian yang and David Zhang, "Metaface learning for +% sparse representation based face recognition", In ICIP, 2010. +%---------------------------------------------------------------------- +% +% Inputs : (1) A : the training data +% (2) X : the coefficient matrix of the training data +% (3) index: the label of the class being processed +% (4) classn: the total number of classes +% (5) Fish_ipts +% . D the dictionary in the last interation +% . trls the labels of training data +% (6) Fish_par +% .dls the labels of the dicitonary atoms +% +% Outputs: (1) Di_new : the updated dictionary of the index-th class +% (2) delet: the indication of deleted atoms +% +%-------------------------------------------------------------------- + +D = Fish_ipts.D; +tau3 = 1; +tau2 = 1; +trls = Fish_ipts.trls; +dls = Fish_par.dls; +delet = []; + +Do = D(:,dls~=index); +Di = D(:,dls==index); +Ai = A(:,trls==index); +Ao = A(:,trls~=index); +Xi = X(:,trls==index); +Xo = X(:,trls~=index); + +Xi_i = Xi(dls==index,:); +Xi_o = Xi(dls~=index,:); +Xo_i = Xo(dls==index,:); +Xo_o = Xo(dls~=index,:); +% X_i = X(dls==index,:); + +Zi = Ai-Do*Xi_o; +Zo = Ao-Do*Xo_o; + +for i = 1: size(Di,2) + Yi = Zi - Di*Xi_i + Di(:,i)*Xi_i(i,:); + Ui = Ai - Di*Xi_i + Di(:,i)*Xi_i(i,:); +% Ua = A - Di*X_i+Di(:,i)*X_i(i,:); + Vo = Zo - Di*Xo_i + Di(:,i)*Xo_i(i,:); + + UaXti = zeros(size(Yi*(Xi_i(i,:))')); + for t_i = 1:classn + if t_i~=index + Xt_i = X(dls==index,trls==t_i); + UaXti = UaXti + (0 - Di*Xt_i + Di(:,i)*Xt_i(i,:))*(Xt_i(i,:))'; + end + end + + tem1 = -Yi*(Xi_i(i,:))' - (tau2)*Ui*(Xi_i(i,:))' - Vo*(Xo_i(i,:))' - tau3*UaXti; +% tem1 = - (tau2)*Ui*(Xi_i(i,:))'; +% tem1 = -Yi*(Xi_i(i,:))' - (tau2)*Ui*(Xi_i(i,:))' - tau3*Ua*(Xo_i(i,:))'; + tem = -tem1; + if norm(tem,2)<1e-6 + Di(:,i) = zeros(size(tem)); + delet = [delet i]; + else + Di(:,i) = tem./norm(tem,2); + end +end + +Di_new = Di; \ No newline at end of file diff --git a/lib/FDDL/utilies/IPM_SC.m b/lib/FDDL/utilies/IPM_SC.m new file mode 100755 index 0000000..fee97a7 --- /dev/null +++ b/lib/FDDL/utilies/IPM_SC.m @@ -0,0 +1,167 @@ +function [opts] = IPM_SC ( ipts , par) +% ======================================================================== +% Spaese Coding of FDDL, Version 1.0 +% Copyright(c) 2011 Meng YANG, Lei Zhang, Xiangchu Feng and David Zhang +% All Rights Reserved. +% +% ----------------------------------------------------------------------- +% Permission to use, copy, or modify this software and its documentation +% for educational and research purposes only and without fee is here +% granted, provided that this copyright notice and the original authors' +% names appear on all copies and supporting documentation. This program +% shall not be used, rewritten, or adapted as the basis of a commercial +% software or hardware product without first obtaining permission of the +% authors. The authors make no representations about the suitability of +% this software for any purpose. It is provided "as is" without express +% or implied warranty. +%---------------------------------------------------------------------- +% +% This is an implementation of the algorithm to do sparse coding +% +% Please refer to the following paper +% +% Meng Yang, Lei Zhang, Xiangchu Feng, and David Zhang,"Fisher Discrimination +% Dictionary Learning for Sparse Representation", In IEEE Int. Conf. on +% Computer Vision, 2011. +% L. Rosasco, A. Verri, M. Santoro, S. Mosci, and S. Villa. Iterative +% Projection Methods for Structured Sparsity Regularization. MIT Technical +% Reports, MIT-CSAIL-TR-2009-050,CBCL-282, 2009. +% J. Bioucas-Dias, M. Figueiredo, ?A new TwIST: two-step iterative shrinkage +% /thresholding algorithms for image restoration?, IEEE Transactions on +% Image Processing, December 2007. +%---------------------------------------------------------------------- +% +% Inputs : (1) ipts : the structre of input data +% .D the dictionary +% .tau1 parameter of l1-norm energy of coefficient +% .y the testing sample +% (2) par : the struture of input parameters +% .eigenv +% +% Outputs: (1) opts : the structure of output data +% .A the coefficient matrix +% .ert the total energy sequence +% +%--------------------------------------------------------------------- + +ipts.cT = 1e+10; % stop criterion +ipts.citeT = 1e-6; % stop criterion +par.nIter = 200; % maximal iteration number +par.cRatio = 1.05; +ipts.initM = 'zero'; % coefficiet initialization method + +switch lower(ipts.initM) + case {'zero'} + x(:,1) = zeros(size(ipts.D,2),1); + case {'transpose'} + x(:,1) = ipts.D'*ipts.y; + case {'pinv'} + x(:,1) = pinv(ipts.D)*ipts.y; + otherwise + error('Nonknown method!'); +end + +A = ipts.D; +tau1 = ipts.tau1; +y = ipts.y; +nIter = par.nIter; +if ~isfield(par,{'eigenv'}) +c = par.cRatio*find_max_eigenv(A'*A); +else +c = par.cRatio*par.eigenv ; +end + +%%%%%%%%%%%%%%%%%%%% +% TWIST parameter +%%%%%%%%%%%%%%%%%%%% +for_ever = 1; +IST_iters = 0; +TwIST_iters = 0; +sparse_sign = 1; +verbose = 1; +enforceMonotone = 1; +lam1 = 1e-4; %default eigenvalues +lamN = 1; %default eigenvalues +rho0 = (1-lam1/lamN)/(1+lam1/lamN); +alpha = 2/(1+sqrt(1-rho0^2)); %default,user can set +beta = alpha*2/(lam1+lamN); %default,user can set + +%%%%%%%%%%%%%%%%%%% +%main loop +%%%%%%%%%%%%%%%%%%% +xm2 = x(:,1); +xm1 = x(:,1); +prev_f = norm(y-A*x(:,1),2)^2+... + tau1*norm(x(:,1),1); +for n_it = 2 : nIter; + + ert(n_it-1) = norm(y-A*x(:,n_it-1),2)^2+... + tau1*norm(x(:,n_it-1),1); + fprintf('Iteration:%f Total gap:%f\n',n_it,ert(n_it-1)); + + while for_ever + % IST(IPM) estimate + v = 1/c*A'*(y-A*xm1)+xm1; + x_temp = soft(v,tau1/c/2); + + if (IST_iters >= 2) | ( TwIST_iters ~= 0) + % set to zero the past when the present is zero + % suitable for sparse inducing priors + if sparse_sign + mask = (x_temp ~= 0); + xm1 = xm1.* mask; + xm2 = xm2.* mask; + end + % two-step iteration + xm2 = (alpha-beta)*xm1 + (1-alpha)*xm2 + beta*x_temp; + % compute residual + f = norm(y-A*xm2,2)^2+... + tau1*norm(xm2,1); + if (f > prev_f) & (enforceMonotone) + TwIST_iters = 0; % do a IST(IPM) iteration if monotonocity fails + else + TwIST_iters = TwIST_iters+1; % TwIST iterations + IST_iters = 0; + x_temp = xm2; + if mod(TwIST_iters,10000) ==0 + c = 0.9*c; + end + break; % break loop while + end + else + f = norm(y-A*x_temp,2)^2+... + tau1*norm(x_temp,1); + if f > prev_f + % if monotonicity fails here is because + % max eig (A'A) > 1. Thus, we increase our guess + % of max_svs + c = 2*c; + if verbose +% fprintf('Incrementing c=%2.2e\n',c); + end + if c > ipts.cT + break; % break loop while + end + IST_iters = 0; + TwIST_iters = 0; + else + TwIST_iters = TwIST_iters + 1; + break; % break loop while + end + end + end + + citerion = abs(f-prev_f)/prev_f; + if citerion < ipts.citeT | c > ipts.cT +% fprintf('Stop!\n c=%2.2e\n citerion=%2.2e\n',c,citerion); + break; + end + + xm2 = xm1; + xm1 = x_temp; + x(:,n_it) = x_temp; + prev_f = f; +end + +opts.x = x(:,end); +opts.ert = ert; \ No newline at end of file diff --git a/lib/FDDL/utilies/soft.m b/lib/FDDL/utilies/soft.m new file mode 100755 index 0000000..d638ce1 --- /dev/null +++ b/lib/FDDL/utilies/soft.m @@ -0,0 +1,4 @@ +% define the soft threshold function, which is used above. +function y = soft(x,tau) + +y = sign(x).*max(abs(x)-tau,0); \ No newline at end of file