naive singular value decomposition draft

svd.rs

    pub fn singular_value_decomp(&self) -> (Matrix<T>, Matrix<T>, Matrix<T>) {
        let ata = self.transpose() * self;
        let (ata_eigenvals, ata_eigvecs) = ata.eigendecomp();

        // (index, √λ) pairs (indices being order of return from `.eigendecomp`)
        let mut enumerated_singular_vals = ata_eigenvals.iter()
            .map(|s| s.sqrt())
            .enumerate()
            .collect::<Vec<_>>();

        enumerated_singular_vals
            // ... in descending order by singular value
            .sort_by(|a, b| b.1.partial_cmp(&a.1).expect("NaN is uncomparable"));

        // Pluck out the eigenvectors in the corresponding order
        let v = enumerated_singular_vals.iter()
            .map(|&(i, _sigma)| ata_eigvecs.select_cols(&[i]))
            .map(|v| { let n = v.norm(); v/n })
            .fold(Matrix::new(ata.rows(), 0, Vec::new()),
                  |vs, vi| vs.hcat(&vi));

        let singular_vals = enumerated_singular_vals.iter()
            .map(|&(_, sigma)| sigma)
            .collect::<Vec<_>>();

        // u₁ = 1/σ₁ Av₁, &c ...
        let u = singular_vals.iter()
            .enumerate()
            // column iterators (issue #52) could make this less awkward
            .map(|(i, sigma)| self * T::one()/sigma * v.select_cols(&[i]))
            .fold(Matrix::new(v.rows(), 0, Vec::new()),
                  |us, ui| us.hcat(&ui));

        (u, Matrix::from_diag(&singular_vals[..]), v.transpose())
    }