summaryrefslogtreecommitdiffstats
path: root/src/armadillo/include/armadillo_bits/op_princomp_meat.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/armadillo/include/armadillo_bits/op_princomp_meat.hpp')
-rw-r--r--src/armadillo/include/armadillo_bits/op_princomp_meat.hpp319
1 files changed, 319 insertions, 0 deletions
diff --git a/src/armadillo/include/armadillo_bits/op_princomp_meat.hpp b/src/armadillo/include/armadillo_bits/op_princomp_meat.hpp
new file mode 100644
index 0000000..db6f83f
--- /dev/null
+++ b/src/armadillo/include/armadillo_bits/op_princomp_meat.hpp
@@ -0,0 +1,319 @@
+// SPDX-License-Identifier: Apache-2.0
+//
+// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
+// Copyright 2008-2016 National ICT Australia (NICTA)
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+// ------------------------------------------------------------------------
+
+
+//! \addtogroup op_princomp
+//! @{
+
+
+
+//! \brief
+//! principal component analysis -- 4 arguments version
+//! computation is done via singular value decomposition
+//! coeff_out -> principal component coefficients
+//! score_out -> projected samples
+//! latent_out -> eigenvalues of principal vectors
+//! tsquared_out -> Hotelling's T^2 statistic
+template<typename T1>
+inline
+bool
+op_princomp::direct_princomp
+ (
+ Mat<typename T1::elem_type>& coeff_out,
+ Mat<typename T1::elem_type>& score_out,
+ Col<typename T1::pod_type>& latent_out,
+ Col<typename T1::elem_type>& tsquared_out,
+ const Base<typename T1::elem_type, T1>& X
+ )
+ {
+ arma_extra_debug_sigprint();
+
+ typedef typename T1::elem_type eT;
+ typedef typename T1::pod_type T;
+
+ const unwrap_check<T1> Y( X.get_ref(), score_out );
+ const Mat<eT>& in = Y.M;
+
+ const uword n_rows = in.n_rows;
+ const uword n_cols = in.n_cols;
+
+ if(n_rows > 1) // more than one sample
+ {
+ // subtract the mean - use score_out as temporary matrix
+ score_out = in; score_out.each_row() -= mean(in);
+
+ // singular value decomposition
+ Mat<eT> U;
+ Col< T> s;
+
+ const bool svd_ok = (n_rows >= n_cols) ? svd_econ(U, s, coeff_out, score_out) : svd(U, s, coeff_out, score_out);
+
+ if(svd_ok == false) { return false; }
+
+ // normalize the eigenvalues
+ s /= std::sqrt( double(n_rows - 1) );
+
+ // project the samples to the principals
+ score_out *= coeff_out;
+
+ if(n_rows <= n_cols) // number of samples is less than their dimensionality
+ {
+ score_out.cols(n_rows-1,n_cols-1).zeros();
+
+ Col<T> s_tmp(n_cols, arma_zeros_indicator());
+
+ s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
+ s = s_tmp;
+
+ // compute the Hotelling's T-squared
+ s_tmp.rows(0,n_rows-2) = T(1) / s_tmp.rows(0,n_rows-2);
+
+ const Mat<eT> S = score_out * diagmat(Col<T>(s_tmp));
+ tsquared_out = sum(S%S,1);
+ }
+ else
+ {
+ // compute the Hotelling's T-squared
+ // TODO: replace with more robust approach
+ const Mat<eT> S = score_out * diagmat(Col<T>( T(1) / s));
+ tsquared_out = sum(S%S,1);
+ }
+
+ // compute the eigenvalues of the principal vectors
+ latent_out = s%s;
+ }
+ else // 0 or 1 samples
+ {
+ coeff_out.eye(n_cols, n_cols);
+
+ score_out.copy_size(in);
+ score_out.zeros();
+
+ latent_out.set_size(n_cols);
+ latent_out.zeros();
+
+ tsquared_out.set_size(n_rows);
+ tsquared_out.zeros();
+ }
+
+ return true;
+ }
+
+
+
+//! \brief
+//! principal component analysis -- 3 arguments version
+//! computation is done via singular value decomposition
+//! coeff_out -> principal component coefficients
+//! score_out -> projected samples
+//! latent_out -> eigenvalues of principal vectors
+template<typename T1>
+inline
+bool
+op_princomp::direct_princomp
+ (
+ Mat<typename T1::elem_type>& coeff_out,
+ Mat<typename T1::elem_type>& score_out,
+ Col<typename T1::pod_type>& latent_out,
+ const Base<typename T1::elem_type, T1>& X
+ )
+ {
+ arma_extra_debug_sigprint();
+
+ typedef typename T1::elem_type eT;
+ typedef typename T1::pod_type T;
+
+ const unwrap_check<T1> Y( X.get_ref(), score_out );
+ const Mat<eT>& in = Y.M;
+
+ const uword n_rows = in.n_rows;
+ const uword n_cols = in.n_cols;
+
+ if(n_rows > 1) // more than one sample
+ {
+ // subtract the mean - use score_out as temporary matrix
+ score_out = in; score_out.each_row() -= mean(in);
+
+ // singular value decomposition
+ Mat<eT> U;
+ Col< T> s;
+
+ const bool svd_ok = (n_rows >= n_cols) ? svd_econ(U, s, coeff_out, score_out) : svd(U, s, coeff_out, score_out);
+
+ if(svd_ok == false) { return false; }
+
+ // normalize the eigenvalues
+ s /= std::sqrt( double(n_rows - 1) );
+
+ // project the samples to the principals
+ score_out *= coeff_out;
+
+ if(n_rows <= n_cols) // number of samples is less than their dimensionality
+ {
+ score_out.cols(n_rows-1,n_cols-1).zeros();
+
+ Col<T> s_tmp(n_cols, arma_zeros_indicator());
+
+ s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2);
+ s = s_tmp;
+ }
+
+ // compute the eigenvalues of the principal vectors
+ latent_out = s%s;
+ }
+ else // 0 or 1 samples
+ {
+ coeff_out.eye(n_cols, n_cols);
+
+ score_out.copy_size(in);
+ score_out.zeros();
+
+ latent_out.set_size(n_cols);
+ latent_out.zeros();
+ }
+
+ return true;
+ }
+
+
+
+//! \brief
+//! principal component analysis -- 2 arguments version
+//! computation is done via singular value decomposition
+//! coeff_out -> principal component coefficients
+//! score_out -> projected samples
+template<typename T1>
+inline
+bool
+op_princomp::direct_princomp
+ (
+ Mat<typename T1::elem_type>& coeff_out,
+ Mat<typename T1::elem_type>& score_out,
+ const Base<typename T1::elem_type, T1>& X
+ )
+ {
+ arma_extra_debug_sigprint();
+
+ typedef typename T1::elem_type eT;
+ typedef typename T1::pod_type T;
+
+ const unwrap_check<T1> Y( X.get_ref(), score_out );
+ const Mat<eT>& in = Y.M;
+
+ const uword n_rows = in.n_rows;
+ const uword n_cols = in.n_cols;
+
+ if(n_rows > 1) // more than one sample
+ {
+ // subtract the mean - use score_out as temporary matrix
+ score_out = in; score_out.each_row() -= mean(in);
+
+ // singular value decomposition
+ Mat<eT> U;
+ Col< T> s;
+
+ const bool svd_ok = (n_rows >= n_cols) ? svd_econ(U, s, coeff_out, score_out) : svd(U, s, coeff_out, score_out);
+
+ if(svd_ok == false) { return false; }
+
+ // project the samples to the principals
+ score_out *= coeff_out;
+
+ if(n_rows <= n_cols) // number of samples is less than their dimensionality
+ {
+ score_out.cols(n_rows-1,n_cols-1).zeros();
+ }
+ }
+ else // 0 or 1 samples
+ {
+ coeff_out.eye(n_cols, n_cols);
+ score_out.copy_size(in);
+ score_out.zeros();
+ }
+
+ return true;
+ }
+
+
+
+//! \brief
+//! principal component analysis -- 1 argument version
+//! computation is done via singular value decomposition
+//! coeff_out -> principal component coefficients
+template<typename T1>
+inline
+bool
+op_princomp::direct_princomp
+ (
+ Mat<typename T1::elem_type>& coeff_out,
+ const Base<typename T1::elem_type, T1>& X
+ )
+ {
+ arma_extra_debug_sigprint();
+
+ typedef typename T1::elem_type eT;
+ typedef typename T1::pod_type T;
+
+ const unwrap<T1> Y( X.get_ref() );
+ const Mat<eT>& in = Y.M;
+
+ if(in.n_elem != 0)
+ {
+ Mat<eT> tmp = in; tmp.each_row() -= mean(in);
+
+ // singular value decomposition
+ Mat<eT> U;
+ Col< T> s;
+
+ const bool svd_ok = (in.n_rows >= in.n_cols) ? svd_econ(U, s, coeff_out, tmp) : svd(U, s, coeff_out, tmp);
+
+ if(svd_ok == false) { return false; }
+ }
+ else
+ {
+ coeff_out.eye(in.n_cols, in.n_cols);
+ }
+
+ return true;
+ }
+
+
+
+template<typename T1>
+inline
+void
+op_princomp::apply
+ (
+ Mat<typename T1::elem_type>& out,
+ const Op<T1,op_princomp>& in
+ )
+ {
+ arma_extra_debug_sigprint();
+
+ const bool status = op_princomp::direct_princomp(out, in.m);
+
+ if(status == false)
+ {
+ out.soft_reset();
+
+ arma_stop_runtime_error("princomp(): decomposition failed");
+ }
+ }
+
+
+
+//! @}
generated by cgit v1.2.1 (git 2.18.0) at 2024-11-22 06:56:00 +0100