pinv broken on column-rank-deficient matrices.

[ixchow]

Applies to mathjs 14.1.0, installed via npm.

May be related to #3012, but that bug does not provide a test case.

According to wikipedia's summary, the pseudo-inverse of a column-rank-deficient matrix mat exists, and acts as a right inverse mat * pinv(mat) = I.

This does not appear to be the case when using mathjs.pinv on matricies with columns that are not linearly independent.

In some cases, the returned matrix is not a right inverse:

const mat = [
	[ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
	[ 0, 0, 1, 1, 0, 0, 0, 0, 0, 1 ],
	[ 0, 0, 2, 0, 1, 0, 0, 1, 0, 0 ],
	[ 1, 0, 1, 0, 0, 0, 0, 0, 1, 0 ],
	[ 0, 1, 1, 0, 1, 1, 2, 0, 0, 1 ]
];

const pinv = math.pinv(mat);
console.table(math.multiply(mat, pinv));

The output here has its first column as all zeros -- quite far from the identity matrix (expected result).

In comparison, computing with wikipedia's summary of how to compute the pseudoinverse works just fine:

const mat = /* same as above */;
//right-inverse recipe, works for matrices with linearly-independent rows:
const pinv = math.multiply(
	math.transpose( mat ),
	math.inv(
		math.multiply(
			mat,
			math.transpose( mat )
		)
	)
);
console.table(math.multiply(mat, pinv));

Here, the output is very close to the identity matrix.

---

Further, pinv just fails on some matrices:

const mat = [
	[ 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
	[ 0, 0, 1, 1, 0, 0, 0, 0, 1 ],
	[ 0, 0, 0, 0, 1, 0, 1, 0, 0 ],
	[ 1, 0, 1, 0, 0, 0, 0, 1, 0 ],
	[ 0, 1, 1, 0, 1, 2, 0, 0, 1 ]
];
const pinv = math.pinv(mat);

Results in an Error: Cannot calculate inverse, determinant is zero exception being thrown. (The right inverse recipe above still works on this matrix.)

---

My suspicion is that the pinv implementation here was only tested on matrices with linearly-independent columns (the rank-surplus case), not linearly-independent rows (the rank-deficient case). It's a useful tool in both situations so it would make sense to patch to deal with both.

[gwhitney]

Thanks! This is a good addition to #3012. Having specific examples is very helpful. You don't happen to have an intuition as how we might find even smaller examples, do you?

[manuelcuesta2005]

Hello, can I work on this issue, I want to know a little more about mathematics in programming.

[gwhitney]

Well, the suggestion in #3012 that this may need a reimplementation of the pseudoinverse in terms of the singular-value decomposition of a matrix. That would involve some reasonably substantial coding, albeit of well-known linear algebra algorithms. You are welcome to give it a shot and submit a PR, but note it will have to pass all of the existing unit tests, and you will need to add unit tests for the above examples as well.

[manuelcuesta2005]

Sounds good to me, so can I start working?

[gwhitney]

I mean, this is a clearly defined bug in mathjs as it stands, so a working PR would be welcome.

[manuelcuesta2005]

Ok, I can understand, but I mean if I have already assigned the issue to work? excuse my insistence or if I am not understanding.

[KeneePatel]

Hi, I’d like to take on this issue. I’ve reviewed the contribution guidelines and done my research on pseudo-inverses, matrix rank, and the concepts of dependent and independent columns. I’ve also looked through the current pinv.js implementation and noticed the TODO about switching to SVD instead of the current rank decomposition. I’m interested in implementing the SVD-based approach, ensuring all tests pass, and updating the documentation accordingly. Please let me know if you’re happy for me to proceed.

[josdejong]

@manuelcuesta2005 thanks for your input. If you want to familiarize yourself with mathematics in programming, it may help to just find one very simple issue in this repository to take on, to get to know the library a bit.

@KeneePatel sounds like you're familiar with the matter which I think is essential, your help is very welcome! Can you pick this issue up? Please let me know if you need any pointer and please keep us up to date if you have any progress (or a lack of it ;) ).

[KeneePatel]

Thank you @josdejong for assigning me with the issue. I will start working upon the issue.

[KeneePatel]

Hello @josdejong,

I've spent the last couple of days digging into SVD and why it's essential for a robust pseudo-inverse implementation. It's clear now that using SVD is necessary for its stability. However, SVD itself is a complex topic. Implementing it efficiently would require quite a bit of dedicated code, and possibly iterative methods.

Given this, I think it makes sense to open a separate issue for SVD first. Once that's in place, implementing a reliable pseudo-inverse should be much more straightforward. What are your thoughts on this?

[gwhitney]

mathjs would welcome a PR simply adding singular value decomposition functionality. It would need to include a thorough suite of unit tests verifying its functionality. If you are so inclined, do feel welcome to submit such a PR.

Thanks very much.

[josdejong]

@KeneePatel please be careful not polluting this issue with AI slop.