2.1.8

README.md

@@ -141,9 +141,6 @@ is a scaled distance, whereby the scaling is derived from the axes of covariance
 * `Cholesky-Banachiewicz matrix decomposition`  
 decomposes any positive definite matrix S (often covariance or comediance matrix) into a product of two triangular matrices: `S = LL'`. The eigenvalues and the determinant are easily obtained from the diagonal of L. We implemented it on `TriangMat` for maximum efficiency. It is used mainly by `eigenvalues`, `eigenvectors`, `mahalanobis` and `pca_reduction`.
 
-* `PCA (Principal Components Analysis)`
-is typically used to reduce the dimensionality of data along some of the most significant directions of variation (eigenvectors).
-
 * `Householder's decomposition`  
 in cases where the precondition (positive definite matrix S) for the Cholesky-Banachiewicz decomposition is not satisfied, Householder's (UR) decomposition is often used as the next best method. It is implemented here on our efficient `struct TriangMat`.
 
@@ -344,7 +341,9 @@ Methods which take an additional generic vector argument, such as a vector of we
 
 ## Appendix: Recent Releases
 
-* **Version 2.1.7** - Removed suspect eigen values/vectors computations. Improved 'cholesky' test.
+* **Version 2.1.8** - Improved `TriangMat::diagonal()`, restored `TriangMat::determinant()`, tidied up `triangmat` test.
+
+* **Version 2.1.7** - Removed suspect eigen values/vectors computations. Improved 'householder' test.
 
 * **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.
 

src/triangmat.rs

@@ -53,19 +53,28 @@ impl TriangMat {
     pub fn diagonal(&self) -> Vec<f64> {
         let mut next = 0_usize;
         let mut skip = 1;
-        self.data
-            .iter()
-            .enumerate()
-            .filter_map(|(i, &x)| {
-                if i == next {
-                    skip += 1;
-                    next = i + skip;
-                    Some(x)
-                } else {
-                    None
-                }
-            })
-            .collect::<Vec<f64>>()
+        let dat = &self.data; 
+        let mut diagonal = Vec::with_capacity(self.dim());
+        while next < dat.len() { 
+            diagonal.push(dat[next]);
+            skip += 1;
+            next += skip;
+        };
+        diagonal
+    }
+    /// Determinant of A = LL' is the square of the product of the diagonal elements
+    /// However, the diagonal elements squared are not the eigenvalues of A!
+    pub fn determinant(&self) -> f64 {
+        let mut next = 0_usize;
+        let mut skip = 1;
+        let mut product = 1_f64;
+        let dat = &self.data; 
+        while next < dat.len() {
+            product *= dat[next];
+            skip += 1;
+            next += skip;  
+        };
+        product*product
     }
     /// New unit (symmetric) TriangMat matrix (data size `n*(n+1)/2`)
     pub fn unit(n: usize) -> Self {
@@ -78,7 +87,7 @@ impl TriangMat {
             data.push(1_f64);
         }
         TriangMat { kind: 2, data }
-    }    
+    }
     /// Translates subscripts to a 1d vector, i.e. natural numbers, to a pair of
     /// (row,column) coordinates within a lower/upper triangular matrix.
     /// Enables memory efficient representation of triangular matrices as one flat vector.

tests/tests.rs

@@ -241,31 +241,29 @@ fn trend() -> Result<(), RE> {
 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 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!("Full Comediance matrix:\n{}",cov.to_full().gr());
+    println!("Symmetric positive definite matrix A:\n{cov}");
+    println!("Full form of A:\n{}",cov.to_full().gr());
     let mut chol = cov.cholesky()?;
-    println!("Cholesky L matrix:\n{chol}"); 
+    println!("Cholesky L matrix, such that A=LL':\n{chol}");
+    println!("Diagonal of L: {}",chol.diagonal().gr());
+    println!("Determinant det(A): {}",chol.determinant().gr());
     let full = chol.to_full();
     chol.transpose();
-    let tfull = chol.to_full();
-    println!("Full L' matrix:\n{}",tfull.gr()); 
-    println!("Comediance reconstructed as LL':\n{}",full.matmult(&tfull)?.gr()); 
+    let tfull = chol.to_full(); 
+    println!("A reconstructed as LL':\n{}",full.matmult(&tfull)?.gr()); 
     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}\
-        \nScale factor: {GR}{:>0.8}{UN}",
+        \nIts Mahalanobis magnitude {GR}{:>8.8}{UN}", 
         dmag,
-        mahamag,
-        mahamag / dmag
+        mahamag
     );
     Ok(())
 }
@@ -288,6 +286,34 @@ fn mat() -> Result<(), RE> {
     println!("\nM'M:\n{}", m.matmult(&t)?.gr());
     Ok(())
 }
+#[test]
+fn householder() -> Result<(), RE> {
+    let a = &[
+        vec![35., 1., 6., 26., 19., 24.],
+        vec![3., 32., 7., 21., 23., 25.],
+        vec![31., 9., 2., 22., 27., 20.],
+        vec![8., 28., 33., 17., 10., 15.],
+        vec![30., 5., 34., 12., 14., 16.],
+        vec![4., 36., 29., 13., 18., 11.],
+    ];
+    let atimesunit = a.matmult(&unit_matrix(a.len()))?;
+    println!("Matrix a:\n{}", atimesunit.gr());
+    let (u, r) = a.house_ur()?;
+    println!("house_ur u' {u}");
+    println!("house_ur r'  {r}");
+    let q = u.house_uapply(&unit_matrix(a.len().min(a[0].len())));
+    println!(
+        "Q matrix\n{}\nOthogonality of Q check (Q'*Q = I):\n{}",
+        q.gr(),
+        q.transpose().matmult(&q)?.gr()
+    );
+    println!("normalized Q;\n{}",q.normalize()?.gr());
+    println!(
+        "Matrix a = QR recreated:\n{}",
+        q.matmult(&r.to_full())?.gr()
+    );
+    Ok(())
+} 
 
 #[test]
 fn vecvec() -> Result<(), RE> {
@@ -520,35 +546,6 @@ fn hulls() -> Result<(), RE> {
     Ok(())
 }
 
-#[test]
-fn householder() -> Result<(), RE> {
-    let a = &[
-        vec![35., 1., 6., 26., 19., 24.],
-        vec![3., 32., 7., 21., 23., 25.],
-        vec![31., 9., 2., 22., 27., 20.],
-        vec![8., 28., 33., 17., 10., 15.],
-        vec![30., 5., 34., 12., 14., 16.],
-        vec![4., 36., 29., 13., 18., 11.],
-    ];
-    let atimesunit = a.matmult(&unit_matrix(a.len()))?;
-    println!("Matrix a:\n{}", atimesunit.gr());
-    let (u, r) = a.house_ur()?;
-    println!("house_ur u' {u}");
-    println!("house_ur r'  {r}");
-    let q = u.house_uapply(&unit_matrix(a.len().min(a[0].len())));
-    println!(
-        "Q matrix\n{}\nOthogonality of Q check (Q'*Q = I):\n{}",
-        q.gr(),
-        q.transpose().matmult(&q)?.gr()
-    );
-    println!("normalized Q;\n{}",q.normalize()?.gr());
-    println!(
-        "Matrix a = QR recreated:\n{}",
-        q.matmult(&r.to_full())?.gr()
-    );
-    Ok(())
-}
-
 #[test]
 fn geometric_medians() -> Result<(), RE> {
     const NAMES: [&str; 5] = [