Add PInv Moore-Penrose pseudoinverse

Open
#292 9 comments 1 reaction 0 assignees View on GitHub

Nobody has claimed this yet.

Assessment

Difficulty
4/5
Estimated time
3-5 days
Newbie friendliness
35/100
Issue type
Feature
Clarity
Mostly clear
Activity status
Stale
Tech stack
rust
Domain
data

Research direction

Start in numrs/src/linalg.rs at the proposed Pinv traits and the existing self.svd(true, true) call. Check how ndarray-linalg represents real and complex Scalar values and how matrix multiplication is supported across those types. Done means a generic Moore-Penrose pseudoinverse with optional rcond that avoids the reported multiplication error and matches the NumPy and MATLAB references.

Written by the indexing model from the issue text.

Description

I was porting some code from NumPy, and needed the pinv[1][2] function (Moore-Penrose pseudoinverse). I have started implementing it but have run into some issues making it generic over real and complex numbers.

pub trait Pinv {
    type B;
    type E;
    type C;
    fn pinv(&self, rcond: Option<Self::E>) -> Result<Self::B, Error>;
}

pub trait PInvVInto {
    type B;
    type E;
    fn pinv(self, rcond: Option<Self::E>) -> Result<Self::B, Error>;
}

pub trait PInvInplace {
    type B;
    type E;
    fn pinv(&mut self, rcond: Option<Self::E>) -> Result<Self::B, Error>;
}

impl<A, S> Pinv for ArrayBase<S, Ix2>
where
    A: Scalar + Lapack + PartialOrd + Zero,
    S: Data<Elem = A>,
{
    type E = A::Real;
    type B = Array2<A>;
    type C = Array1<A::Real>;

    fn pinv(&self, rconf: Option<Self::E>) -> Result<Self::B, Error> {
        let result: (Option<Self::B>, Self::C, Option<Self::B>) =
            self.svd(true, true).map_err(|_| Error::SvdFailed)?;
        if let (Some(u), s, Some(vt)) = result {
            // discard small singular values
            let s: Self::C = if let Some(rconf) = rconf {
                s.mapv(|v| if v > rconf { v } else { Zero::zero() })
            } else {
                s
            };

            let u = u.reversed_axes();
            let s = s.insert_axis(Axis(1));
            let x1 = s * u;
            Ok(&vt.reversed_axes() * x1)
        } else {
            Err(Error::SvdFailed)
        }
    }
}

At the moment I am getting this error

error[E0277]: cannot multiply `ndarray::ArrayBase<OwnedRepr<<A as ndarray_linalg::Scalar>::Real>, Dim<[usize; 2]>>` by `ndarray::ArrayBase<OwnedRepr<A>, Dim<[usize; 2]>>`
  --> numrs/src/linalg.rs:58:24
   |
58 |             let x1 = s * u;
   |                        ^ no implementation for `ndarray::ArrayBase<OwnedRepr<<A as ndarray_linalg::Scalar>::Real>, Dim<[usize; 2]>> * ndarray::ArrayBase<OwnedRepr<A>, Dim<[usize; 2]>>`
   |
   = help: the trait `Mul<ndarray::ArrayBase<OwnedRepr<A>, Dim<[usize; 2]>>>` is not implemented for `ndarray::ArrayBase<OwnedRepr<<A as ndarray_linalg::Scalar>::Real>, Dim<[usize; 2]>>`
help: consider extending the `where` bound, but there might be an alternative better way to express this requirement
   |
38 |     S: Data<Elem = A>, ndarray::ArrayBase<OwnedRepr<<A as ndarray_linalg::Scalar>::Real>, Dim<[usize; 2]>>: Mul<ndarray::ArrayBase<OwnedRepr<A>, Dim<[usize; 2]>>>

My understanding this is because I am trying to multiple Scalar::Real with a Scalar, but I am not sure how this should be fixed.

  1. https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html
  2. https://www.mathworks.com/help/matlab/ref/pinv.html
Dominant language
Rust
Stars
452
Forks
95
PR merge metrics
No merged PRs in 30d

Contributor guide

No contributing guide indexed for this repository

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from rust-ndarray/ndarray-linalg

All issues in rust-ndarray/ndarray-linalg

Similar issues

More Rust issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.