2.1.6

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "rstats"
-version = "2.1.5"
+version = "2.1.6"
 authors = ["Libor Spacek"]
 edition = "2021"
 description = "Statistics, Information Measures, Data Analysis, Linear Algebra, Clifford Algebra, Machine Learning, Geometric Median, Matrix Decompositions, PCA, Mahalanobis Distance, Hulls, Multithreading.."

README.md

@@ -344,11 +344,13 @@ Methods which take an additional generic vector argument, such as a vector of we
 
 ## Appendix: Recent Releases
 
-* **Version 2.1.5** - Added `projection` to trait `VecVecg` to project all self vectors to a new basis. This can be used e.g. for Principal Components Analysis data reduction, using some chosen eigenvectors as the new basis.
+* **Version 2.1.6** - Added `eigen` and `principals` to TriangMat. `Eigen` computes eigenvalues and eigenvectors using Householder's UR decomposition. `Principals` uses these to order the eigenvectors by absolute values of eigenvalues and to choose the requested number of principal components.
+
+* **Version 2.1.5** - Added `projection` to trait `VecVecg` to project all self vectors to a new basis. This can be used e.g. for Principal Components Analysis data reduction, using some of the eigenvectors as the new basis.
 
 * **Version 2.1.4** - Tidied up some error processing.
 
-* **Version 2.1.3** - Changed `eigenvectors` to compute normalized eigenvectors of the original data rather than of its covariance matrix. That is now done by better named `normalize` (should you still need it). `Eigenvectors` solves forward substitution to find each vector.
+* **Version 2.1.3** - Added `normalize` (normalize columns of a matrix and transpose them to rows).
 
 * **Version 2.1.2** - Added function `project` to project a `TriangMat` to a lower dimensional space of selected dimensions. Removed `rows` which was a duplicate of `dim`.
 
@@ -358,8 +360,6 @@ Methods which take an additional generic vector argument, such as a vector of we
 
 * **Version 2.0.11** - removed not so useful `variances`. Tidied up error processing in `vecvecg.rs`. Added to it `serial_covar` and `serial_wcovar` for when heavy loading of all the cores may not be wanted.
 
-* **Version 2.0.10** - Added `eigenvectors` to `struct TriangMat` (to enable PCA).
-
 * **Version 2.0.9** - Pruned some rarely used methods, simplified `gmparts` and `gmerror`, updated dependencies.
 
 * **Version 2.0.8**' - Changed initial guess in iterative weighted gm methods to weighted mean. This, being more accurate than plain mean, leads to fewer iterations. Updated some dependencies.

src/lib.rs

@@ -289,10 +289,12 @@ pub trait VecVec<T> {
     fn column(self, cnum: usize) -> Vec<f64>;
     /// Transpose slice of vecs matrix
     fn transpose(self) -> Vec<Vec<f64>>;
-    /// Normalize columns, so that they are all unit vectors
+    /// Normalize columns of a matrix and transpose them to rows
     fn normalize(self) -> Result<Vec<Vec<f64>>, RE>;
+    /// Projects self onto a given basis, e.g. PCA dimensional reduction. 
+    fn projection(self, basis: &[Vec<f64>]) -> Result<Vec<Vec<f64>>, RE>;
     /// Householder's method returning matrices (U,R)
-    fn house_ur(self) -> Result<(TriangMat, TriangMat), RE>;
+    fn house_ur(self) -> Result<(TriangMat, TriangMat), RE>; 
     /// Joint probability density function of n matched slices of the same length
     fn jointpdfn(self) -> Result<Vec<f64>, RE>;
     /// Joint entropy between a set of vectors of the same length
@@ -347,12 +349,10 @@ pub trait VecVec<T> {
     fn par_gmedian(self, eps: f64) -> Vec<f64>;
     /// Like `gmedian` but returns also the sum of reciprocals of distances
     fn gmparts(self, eps: f64) -> (Vec<f64>, f64);
-    /// Flattened lower triangular part of a covariance matrix of a Vec of f64 vectors.
+    /// Lower triangular part of a (symmetric) covariance matrix of a Vec of f64 vectors.
     fn covar(self, mid: &[f64]) -> Result<TriangMat, RE>;
-    /// Flattened lower triangular part of a covariance matrix of a Vec of f64 vectors.
-    fn serial_covar(self, mid: &[f64]) -> Result<TriangMat, RE>;
-    /// Projects self onto a given basis, e.g. PCA dimensional reduction. 
-    fn projection(self, basis: &[Vec<f64>]) -> Result<Vec<Vec<f64>>, RE>;
+    /// Lower triangular part of a (symmetric) covariance matrix of a Vec of f64 vectors.
+    fn serial_covar(self, mid: &[f64]) -> Result<TriangMat, RE>; 
 }
 
 /// Methods applicable to slice of vectors of generic end type, plus one other argument

src/triangmat.rs

@@ -79,50 +79,31 @@ impl TriangMat {
         }
         TriangMat { kind: 2, data }
     }
-    /// Eigenvalues from Cholesky L matrix
-    pub fn eigenvalues(&self) -> Vec<f64> {
-        self.diagonal().iter().map(|&x| x * x).collect::<Vec<f64>>()
+
+    /// Eigenvalues and normalised eigenvectors of a symmetric (covariance matrix) C,
+    /// using Householder's decomposition C=UR.
+    pub fn eigen(&self) ->  Result<(Vec<f64>,Vec<Vec<f64>>), RE> {
+        let (u, r) = self.to_full().house_ur()?; 
+        Ok((
+            r.diagonal().iter().map(|&e| e.abs()).collect(),
+            u.house_uapply(&unit_matrix(r.dim())).normalize()?))
     }
-    /// Determinant from Cholesky L matrix
-    pub fn determinant(&self) -> f64 {
-        self.diagonal().iter().map(|&x| x * x).product()
-    }
-    /// Normalized full rows from a triangular matrix.
-    /// When the matrix is symmetric, e.g. a covariance matrix, the result are its normalized eigenvectors
-    pub fn normalize(&self) -> Vec<Vec<f64>> {
-        let mut fullcov = self.to_full();
-        fullcov
-            .iter_mut()
-            .for_each(|eigenvector| eigenvector.munit());
-        fullcov
-    }
-    /// Normalized eigenvectors of A=LL'.
-    /// Where L, supplied as self, is the lower triangular Cholesky decomposition 
-    /// of covariance/comediance matrix for A.
+
+    /// Principal components of symmetric lower triangular matrix, such as covariance
     /// Returns `choose` number of eigenvectors corresponding to the largest eigenvalues.
-    pub fn eigenvectors(&self, choose: usize) -> Result<Vec<Vec<f64>>, RE> {
+    pub fn principals(&self, choose: usize) -> Result<Vec<Vec<f64>>, RE> {
         if self.is_empty() {
             return nodata_error("eigenvectors: empty L");
         };
         let n = self.dim();
         if choose > n {
             return data_error("eigenvectors: choose is more than the number of eigenvectors");
         };
-        let mut eigenvectors = Vec::new();
-        let eigenvals = self.eigenvalues();
-        for (rownum, &eigenvalue) in eigenvals.iter().enumerate() {
-            let mut padded_row = self.row(rownum);
-            for _i in rownum+1..n {
-                padded_row.push(0_f64);
-            }
-            let mut eigenvec = self.forward_substitute(&padded_row.smult(eigenvalue))?;
-            eigenvec.munit(); // normalize the eigenvector
-            eigenvectors.push(eigenvec);
-        }
+        let (eigenvals,eigenvecs) = self.eigen()?; 
         // descending sort index of eigenvalues
         let mut index = eigenvals.isort_indexed(0..n, |a, b| b.total_cmp(a));
         index.truncate(choose); // keep only `choose` best
-        let pruned = index.unindex(&eigenvectors, true);
+        let pruned = index.unindex(&eigenvecs, true);
         Ok(pruned)
     }
 

src/vecvec.rs

@@ -57,7 +57,7 @@ where
             .collect()
     }
 
-    /// Normalize columns, so that they become unit row vectors
+    /// Normalize columns of a matrix, so that they become unit row vectors
     fn normalize(self) -> Result<Vec<Vec<f64>>, RE> {
         (0..self[0].len())
             .into_par_iter()
@@ -86,7 +86,7 @@ where
                 rlast.extend(rvec);
                 // drained, reflected with this uvec, and rebuilt, all remaining rows of r
             }
-            // these uvecs are columns, so they must saved column-wise
+            // these uvecs are columns, so they must be saved column-wise
             for (row, &usave) in uvec.iter().enumerate() {
                 ures[sumn(row + j) + j] = usave; // using triangular index
             }

tests/tests.rs

@@ -242,24 +242,30 @@ fn triangmat() -> Result<(), RE> {
     println!("\n{}", TriangMat::unit(7).gr());
     println!("\n{}", TriangMat::unit(7).to_full().gr());   
     println!("Diagonal: {}",TriangMat::unit(7).diagonal().gr());
-    let d = 10_usize;
-    let n = 90_usize;
+    // let cov = pts.covar(&pts.par_gmedian(EPS))?;
+    let cov = TriangMat{kind:2,
+        data:vec![2_f64,1.,3.,0.5,0.2,1.5,0.3,0.1,0.4,2.5]};
+    println!("Comediance matrix:\n{cov}");
+    println!("Projected to subspace given by [0,2,3]:\n{}",cov.project(&[0,2,3]));
+    let (evs,evecs) = cov.eigen()?;
+    println!("Eigenvalues by Householder decomposition:\n{}",evs.gr());
+    println!("Determinant (their product): {}",evs.iter().product::<f64>().gr()); 
+    println!("Eigenvectors by Householder decomposition:\n{}",evecs.gr());
+    let principals = cov.principals(2)?;
+    println!("Two Principal Components:\n{}",principals.gr());
+    let n = 10_usize;
+    let d = 4;
     println!("Testing on a random set of {n} points in {d} dimensional space"); 
     let pts = ranvv_f64_range(n,d, 0.0..=4.0)?;
-    // println!("\nTest data:\n{}",pts.gr());
-    // let transppt = pts.transpose();
-    let cov = pts.covar(&pts.par_gmedian(EPS))?;
-    println!("Comediance matrix:\n{cov}");
-    println!("Projected to subspace given by [0,2,4,6,9]:\n{}",cov.project(&[0,2,4,6,9]));
+    let newdat = pts.projection(&principals)?; 
+    println!("Random data projected (PCA reduced) onto these two eigenvectors:\n{}\n",newdat.gr());
     let chol = cov.cholesky()?;
-    println!("Cholesky L matrix:\n{chol}");
-    println!("Eigenvalues by Cholesky decomposition:\n{}",
-        chol.eigenvalues().gr());
-    println!("Determinant (their product): {}",chol.determinant().gr());  
-    let d = ranv_f64(d)?;
-    let dmag = d.vmag();
-    let mahamag = chol.mahalanobis(&d)?;
-    println!("Random test vector:\n{}", d.gr());
+    println!("Cholesky L matrix:\n{chol}"); 
+    let d = 4_usize;
+    let dv = ranv_f64(d)?;
+    let dmag = dv.vmag();
+    let mahamag = chol.mahalanobis(&dv)?;
+    println!("Random test vector:\n{}", dv.gr());
     println!(
         "Its Euclidian magnitude   {GR}{:>8.8}{UN}\
         \nIts Mahalanobis magnitude {GR}{:>8.8}{UN}\
@@ -268,10 +274,6 @@ fn triangmat() -> Result<(), RE> {
         mahamag,
         mahamag / dmag
     );
-    let evecs = chol.eigenvectors(3)?;
-    println!("Best three unit eigenvectors:\n{}",evecs.gr());
-    let newdat = pts.projection(&evecs)?; 
-    println!("Original data projected (PCA reduced) onto these eigenvectors:\n{}",newdat.gr());
     Ok(())
 }