op_princomp_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Dimitrios Bouzas (dimitris dot mpouzas at gmail dot com)
00006 // - Conrad Sanderson (conradsand at ieee dot org)
00007 // 
00008 // This file is part of the Armadillo C++ library.
00009 // It is provided without any warranty of fitness
00010 // for any purpose. You can redistribute this file
00011 // and/or modify it under the terms of the GNU
00012 // Lesser General Public License (LGPL) as published
00013 // by the Free Software Foundation, either version 3
00014 // of the License or (at your option) any later version.
00015 // (see http://www.opensource.org/licenses for more info)
00016 
00017 
00018 //! \addtogroup op_princomp
00019 //! @{
00020 
00021 
00022 
00023 //! \brief
00024 //! principal component analysis -- 4 arguments version
00025 //! computation is done via singular value decomposition
00026 //! coeff_out    -> principal component coefficients
00027 //! score_out    -> projected samples
00028 //! latent_out   -> eigenvalues of principal vectors
00029 //! tsquared_out -> Hotelling's T^2 statistic
00030 template<typename eT>
00031 inline
00032 void
00033 op_princomp::direct_princomp
00034   (
00035         Mat<eT>& coeff_out,
00036         Mat<eT>& score_out,
00037         Col<eT>& latent_out, 
00038         Col<eT>& tsquared_out,
00039   const Mat<eT>& in
00040   )
00041   {
00042   arma_extra_debug_sigprint();
00043 
00044   const u32 n_rows = in.n_rows;
00045   const u32 n_cols = in.n_cols;
00046   
00047   if(n_rows > 1) // more than one sample
00048     {
00049     // subtract the mean - use score_out as temporary matrix
00050     score_out = in - repmat(mean(in), n_rows, 1);
00051     
00052     // singular value decomposition
00053     Mat<eT> U;
00054     Col<eT> s;
00055     
00056     const bool svd_ok = svd(U,s,coeff_out,score_out);
00057     
00058     if(svd_ok == false)
00059       {
00060       arma_print("princomp(): singular value decomposition failed");
00061       
00062       coeff_out.reset();
00063       score_out.reset();
00064       latent_out.reset();
00065       tsquared_out.reset();
00066       
00067       return;
00068       }
00069     
00070     
00071     //U.reset();  // TODO: do we need this ?  U will get automatically deleted anyway
00072     
00073     // normalize the eigenvalues
00074     s /= std::sqrt(n_rows - 1);
00075     
00076     // project the samples to the principals
00077     score_out *= coeff_out;
00078     
00079     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00080       {
00081       score_out.cols(n_rows-1,n_cols-1).zeros();
00082       
00083       //Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
00084       Col<eT> s_tmp(n_cols);
00085       s_tmp.zeros();
00086       
00087       s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00088       s = s_tmp;
00089           
00090       // compute the Hotelling's T-squared
00091       s_tmp.rows(0,n_rows-2) = eT(1) / s_tmp.rows(0,n_rows-2);
00092       
00093       const Mat<eT> S = score_out * diagmat(Col<eT>(s_tmp));   
00094       tsquared_out = sum(S%S,1); 
00095       }
00096     else
00097       {
00098       // compute the Hotelling's T-squared   
00099       const Mat<eT> S = score_out * diagmat(Col<eT>( eT(1) / s));
00100       tsquared_out = sum(S%S,1);
00101       }
00102             
00103     // compute the eigenvalues of the principal vectors
00104     latent_out = s%s;
00105     }
00106   else // single sample - row
00107     {
00108     if(n_rows == 1)
00109       {
00110       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00111       
00112       score_out.copy_size(in);
00113       score_out.zeros();
00114       
00115       latent_out.set_size(n_cols);
00116       latent_out.zeros();
00117       
00118       tsquared_out.set_size(1);
00119       tsquared_out.zeros();    
00120       }
00121     else
00122       {
00123       coeff_out.reset();
00124       score_out.reset();
00125       latent_out.reset();
00126       tsquared_out.reset();
00127       }
00128     }
00129   
00130   }
00131 
00132 
00133 
00134 //! \brief
00135 //! principal component analysis -- 3 arguments version
00136 //! computation is done via singular value decomposition
00137 //! coeff_out    -> principal component coefficients
00138 //! score_out    -> projected samples
00139 //! latent_out   -> eigenvalues of principal vectors
00140 template<typename eT>
00141 inline
00142 void
00143 op_princomp::direct_princomp
00144   (
00145         Mat<eT>& coeff_out,
00146         Mat<eT>& score_out,
00147         Col<eT>& latent_out,
00148   const Mat<eT>& in
00149   )
00150   {
00151   arma_extra_debug_sigprint();
00152   
00153   const u32 n_rows = in.n_rows;
00154   const u32 n_cols = in.n_cols;
00155   
00156   if(n_rows > 1) // more than one sample
00157     {
00158     // subtract the mean - use score_out as temporary matrix
00159     score_out = in - repmat(mean(in), n_rows, 1);
00160     
00161     // singular value decomposition
00162     Mat<eT> U;
00163     Col<eT> s;
00164     
00165     const bool svd_ok = svd(U,s,coeff_out,score_out);
00166     
00167     if(svd_ok == false)
00168       {
00169       arma_print("princomp(): singular value decomposition failed");
00170       
00171       coeff_out.reset();
00172       score_out.reset();
00173       latent_out.reset();
00174       
00175       return;
00176       }
00177     
00178     
00179     // U.reset();
00180     
00181     // normalize the eigenvalues
00182     s /= std::sqrt(n_rows - 1);
00183     
00184     // project the samples to the principals
00185     score_out *= coeff_out;
00186     
00187     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00188       {
00189       score_out.cols(n_rows-1,n_cols-1).zeros();
00190       
00191       Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
00192       s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00193       s = s_tmp;
00194       }
00195       
00196     // compute the eigenvalues of the principal vectors
00197     latent_out = s%s;
00198     
00199     }
00200   else // single sample - row
00201     {
00202     if(n_rows == 1)
00203       {
00204       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00205       
00206       score_out.copy_size(in);
00207       score_out.zeros();
00208       
00209       latent_out.set_size(n_cols);
00210       latent_out.zeros(); 
00211       }
00212     else
00213       {
00214       coeff_out.reset();
00215       score_out.reset();
00216       latent_out.reset();
00217       }
00218     }
00219   
00220   }
00221 
00222 
00223 
00224 //! \brief
00225 //! principal component analysis -- 2 arguments version
00226 //! computation is done via singular value decomposition
00227 //! coeff_out    -> principal component coefficients
00228 //! score_out    -> projected samples
00229 template<typename eT>
00230 inline
00231 void
00232 op_princomp::direct_princomp
00233   (
00234         Mat<eT>& coeff_out,
00235         Mat<eT>& score_out,
00236   const Mat<eT>& in
00237   )
00238   {
00239   arma_extra_debug_sigprint();
00240 
00241   const u32 n_rows = in.n_rows;
00242   const u32 n_cols = in.n_cols;
00243   
00244   if(n_rows > 1) // more than one sample
00245     {
00246     // subtract the mean - use score_out as temporary matrix
00247     score_out = in - repmat(mean(in), n_rows, 1);
00248     
00249     // singular value decomposition
00250     Mat<eT> U;
00251     Col<eT> s;
00252     
00253     const bool svd_ok = svd(U,s,coeff_out,score_out);
00254     
00255     if(svd_ok == false)
00256       {
00257       arma_print("princomp(): singular value decomposition failed");
00258       
00259       coeff_out.reset();
00260       score_out.reset();
00261       
00262       return;
00263       }
00264     
00265     // U.reset();
00266     
00267     // normalize the eigenvalues
00268     s /= std::sqrt(n_rows - 1);
00269     
00270     // project the samples to the principals
00271     score_out *= coeff_out;
00272     
00273     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00274       {
00275       score_out.cols(n_rows-1,n_cols-1).zeros();
00276       
00277       Col<eT> s_tmp = zeros< Col<eT> >(n_cols);
00278       s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00279       s = s_tmp;
00280       }
00281     }
00282   else // single sample - row
00283     {
00284     if(n_rows == 1)
00285       {
00286       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00287       score_out.copy_size(in);
00288       score_out.zeros();
00289       }
00290     else
00291       {
00292       coeff_out.reset();
00293       score_out.reset();
00294       }
00295     }
00296   }
00297 
00298 
00299 
00300 //! \brief
00301 //! principal component analysis -- 1 argument version
00302 //! computation is done via singular value decomposition
00303 //! coeff_out    -> principal component coefficients
00304 template<typename eT>
00305 inline
00306 void
00307 op_princomp::direct_princomp
00308   (
00309         Mat<eT>& coeff_out,
00310   const Mat<eT>& in
00311   )
00312   {
00313   arma_extra_debug_sigprint();
00314   
00315   if(in.n_elem != 0)
00316     {
00317     // singular value decomposition
00318     Mat<eT> U;
00319     Col<eT> s;
00320     
00321     const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
00322     
00323     const bool svd_ok = svd(U,s,coeff_out, tmp);
00324     
00325     if(svd_ok == false)
00326       {
00327       arma_print("princomp(): singular value decomposition failed");
00328       
00329       coeff_out.reset();
00330       }
00331     }
00332   else
00333     {
00334     coeff_out.reset();
00335     }
00336   }
00337 
00338 
00339 
00340 //! \brief
00341 //! principal component analysis -- 4 arguments complex version
00342 //! computation is done via singular value decomposition
00343 //! coeff_out    -> principal component coefficients
00344 //! score_out    -> projected samples
00345 //! latent_out   -> eigenvalues of principal vectors
00346 //! tsquared_out -> Hotelling's T^2 statistic
00347 template<typename T>
00348 inline
00349 void
00350 op_princomp::direct_princomp
00351   (
00352         Mat< std::complex<T> >& coeff_out,
00353         Mat< std::complex<T> >& score_out,
00354         Col<T>&                 latent_out,
00355         Col< std::complex<T> >& tsquared_out,
00356   const Mat< std::complex<T> >& in
00357   )
00358   {
00359   arma_extra_debug_sigprint();
00360   
00361   typedef std::complex<T> eT;
00362   
00363   const u32 n_rows = in.n_rows;
00364   const u32 n_cols = in.n_cols;
00365   
00366   if(n_rows > 1) // more than one sample
00367     {
00368     // subtract the mean - use score_out as temporary matrix
00369     score_out = in - repmat(mean(in), n_rows, 1);
00370     
00371     // singular value decomposition
00372     Mat<eT> U;
00373     Col<T> s;
00374     
00375     const bool svd_ok = svd(U,s,coeff_out,score_out); 
00376     
00377     if(svd_ok == false)
00378       {
00379       arma_print("princomp(): singular value decomposition failed");
00380       
00381       coeff_out.reset();
00382       score_out.reset();
00383       latent_out.reset();
00384       tsquared_out.reset();
00385       
00386       return;
00387       }
00388     
00389     
00390     //U.reset();
00391     
00392     // normalize the eigenvalues
00393     s /= std::sqrt(n_rows - 1);
00394     
00395     // project the samples to the principals
00396     score_out *= coeff_out;
00397     
00398     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00399       {
00400       score_out.cols(n_rows-1,n_cols-1).zeros();
00401       
00402       Col<T> s_tmp = zeros< Col<T> >(n_cols);
00403       s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00404       s = s_tmp;
00405           
00406       // compute the Hotelling's T-squared   
00407       s_tmp.rows(0,n_rows-2) = 1.0 / s_tmp.rows(0,n_rows-2);
00408       const Mat<eT> S = score_out * diagmat(Col<T>(s_tmp));                     
00409       tsquared_out = sum(S%S,1); 
00410       }
00411     else
00412       {
00413       // compute the Hotelling's T-squared   
00414       const Mat<eT> S = score_out * diagmat(Col<T>(T(1) / s));                     
00415       tsquared_out = sum(S%S,1);
00416       }
00417     
00418     // compute the eigenvalues of the principal vectors
00419     latent_out = s%s;
00420     
00421     }
00422   else // single sample - row
00423     {
00424     if(n_rows == 1)
00425       {
00426       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00427       
00428       score_out.copy_size(in);
00429       score_out.zeros();
00430       
00431       latent_out.set_size(n_cols);
00432       latent_out.zeros();
00433       
00434       tsquared_out.set_size(1);
00435       tsquared_out.zeros();    
00436       }
00437     else
00438       {
00439       coeff_out.reset();
00440       score_out.reset();
00441       latent_out.reset();
00442       tsquared_out.reset();
00443       }
00444     }
00445   }
00446 
00447 
00448 
00449 //! \brief
00450 //! principal component analysis -- 3 arguments complex version
00451 //! computation is done via singular value decomposition
00452 //! coeff_out    -> principal component coefficients
00453 //! score_out    -> projected samples
00454 //! latent_out   -> eigenvalues of principal vectors
00455 template<typename T>
00456 inline
00457 void
00458 op_princomp::direct_princomp
00459   (
00460         Mat< std::complex<T> >& coeff_out,
00461         Mat< std::complex<T> >& score_out,
00462         Col<T>&                 latent_out,
00463   const Mat< std::complex<T> >& in
00464   )
00465   {
00466   arma_extra_debug_sigprint();
00467   
00468   typedef std::complex<T> eT;
00469   
00470   const u32 n_rows = in.n_rows;
00471   const u32 n_cols = in.n_cols;
00472   
00473   if(n_rows > 1) // more than one sample
00474     {
00475     // subtract the mean - use score_out as temporary matrix
00476     score_out = in - repmat(mean(in), n_rows, 1);
00477     
00478     // singular value decomposition
00479     Mat<eT> U;
00480     Col< T> s;
00481     
00482     const bool svd_ok = svd(U,s,coeff_out,score_out);
00483     
00484     if(svd_ok == false)
00485       {
00486       arma_print("princomp(): singular value decomposition failed");
00487       
00488       coeff_out.reset();
00489       score_out.reset();
00490       latent_out.reset();
00491       
00492       return;
00493       }
00494     
00495     
00496     // U.reset();
00497     
00498     // normalize the eigenvalues
00499     s /= std::sqrt(n_rows - 1);
00500     
00501     // project the samples to the principals
00502     score_out *= coeff_out;
00503     
00504     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00505       {
00506       score_out.cols(n_rows-1,n_cols-1).zeros();
00507       
00508       Col<T> s_tmp = zeros< Col<T> >(n_cols);
00509       s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
00510       s = s_tmp;
00511       }
00512       
00513     // compute the eigenvalues of the principal vectors
00514     latent_out = s%s;
00515 
00516     }
00517   else // single sample - row
00518     {
00519     if(n_rows == 1)
00520       {
00521       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00522       score_out.copy_size(in);
00523       score_out.zeros();
00524       latent_out.set_size(n_cols);
00525       latent_out.zeros();
00526       }
00527     else
00528       {
00529       coeff_out.reset();
00530       score_out.reset();
00531       latent_out.reset();
00532       }
00533     }
00534   }
00535 
00536 
00537 
00538 //! \brief
00539 //! principal component analysis -- 2 arguments complex version
00540 //! computation is done via singular value decomposition
00541 //! coeff_out    -> principal component coefficients
00542 //! score_out    -> projected samples
00543 template<typename T>
00544 inline
00545 void
00546 op_princomp::direct_princomp
00547   (
00548         Mat< std::complex<T> >& coeff_out,
00549         Mat< std::complex<T> >& score_out,
00550   const Mat< std::complex<T> >& in
00551   )
00552   {
00553   arma_extra_debug_sigprint();
00554   
00555   typedef std::complex<T> eT;
00556   
00557   const u32 n_rows = in.n_rows;
00558   const u32 n_cols = in.n_cols;
00559   
00560   if(n_rows > 1) // more than one sample
00561     {
00562     // subtract the mean - use score_out as temporary matrix
00563     score_out = in - repmat(mean(in), n_rows, 1);
00564     
00565     // singular value decomposition
00566     Mat<eT> U;
00567     Col< T> s;
00568     
00569     const bool svd_ok = svd(U,s,coeff_out,score_out);
00570     
00571     if(svd_ok == false)
00572       {
00573       arma_print("princomp(): singular value decomposition failed");
00574       
00575       coeff_out.reset();
00576       score_out.reset();
00577       
00578       return;
00579       }
00580     
00581     // U.reset();
00582     
00583     // normalize the eigenvalues
00584     s /= std::sqrt(n_rows - 1);
00585 
00586     // project the samples to the principals
00587     score_out *= coeff_out;
00588 
00589     if(n_rows <= n_cols) // number of samples is less than their dimensionality
00590       {
00591       score_out.cols(n_rows-1,n_cols-1).zeros();
00592       }
00593 
00594     }
00595   else // single sample - row
00596     {
00597     if(n_rows == 1)
00598       {
00599       coeff_out = eye< Mat<eT> >(n_cols, n_cols);
00600       
00601       score_out.copy_size(in);
00602       score_out.zeros();
00603       }
00604     else
00605       {
00606       coeff_out.reset();
00607       score_out.reset();
00608       }
00609     }
00610   }
00611 
00612 
00613 
00614 //! \brief
00615 //! principal component analysis -- 1 argument complex version
00616 //! computation is done via singular value decomposition
00617 //! coeff_out    -> principal component coefficients
00618 template<typename T>
00619 inline
00620 void
00621 op_princomp::direct_princomp
00622   (
00623         Mat< std::complex<T> >& coeff_out,
00624   const Mat< std::complex<T> >& in
00625   )
00626   {
00627   arma_extra_debug_sigprint();
00628   
00629   typedef typename std::complex<T> eT;
00630   
00631   if(in.n_elem != 0)
00632     {
00633     // singular value decomposition
00634     Mat<eT> U;
00635     Col< T> s;
00636     
00637     const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1);
00638     
00639     const bool svd_ok = svd(U,s,coeff_out, tmp);
00640     
00641     if(svd_ok == false)
00642       {
00643       arma_print("princomp(): singular value decomposition failed");
00644       
00645       coeff_out.reset();
00646       }
00647     }
00648   else
00649     {
00650     coeff_out.reset();
00651     }
00652   }
00653 
00654 
00655 
00656 template<typename T1>
00657 inline
00658 void
00659 op_princomp::apply
00660   (
00661         Mat<typename T1::elem_type>& out,
00662   const Op<T1,op_princomp>&          in
00663   )
00664   {
00665   arma_extra_debug_sigprint();
00666   
00667   typedef typename T1::elem_type eT;
00668   
00669   const unwrap_check<T1> tmp(in.m, out);
00670   const Mat<eT>& A     = tmp.M;
00671   
00672   op_princomp::direct_princomp(out, A);
00673   }
00674 
00675 
00676 
00677 //! @}