Matrix inversion in Ivy

[fzipp]

Should anyone require it, here's an implementation of matrix inversion via Gauss-Jordan elimination in Ivy, similar to this APL implementation.

# Identity matrix
op id n = (iota n) o.== iota n

op i elim M =
	a = abs i drop M[; i]
	p = i + first a iota max/ a
	M[i p] = M[p i]
	M[i] = M[i] / M[i; i]
	M - (M[; i] * i != iota count M) o.* M[i; ]

# Matrix division
op A mdiv B =
	(rho A) rho (0 1 * rho B) drop elim/ (rot iota min/ rho B) , box B , A

# Matrix inverse
op inv M =
	(id count M) mdiv M

Example usage:

M = 4 4 rho 0 1 2 -2 0 2 1 -1 0 2 -2 1 1 0 5 4; M
 0  1  2 -2
 0  2  1 -1
 0  2 -2  1
 1  0  5  4

inv M
 44/3 -49/3     9     1
 -1/3   2/3     0     0
 -4/3   5/3    -1     0
   -2     2    -1     0

M +.* inv M
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

[robpike]

Thanks. If someone didn't keep asking for other features :), my plan was to implement it this week.

But: I want to do Moore-Penrose general inversion, although Gauss-Jordan elimination is one of the steps for that.

By the way, are you just playing with ivy or do you have a real use for it?

[fzipp]

A built-in implementation would, of course, be even better.

By the way, are you just playing with ivy or do you have a real use for it?

Mainly, I play around with it, learn, and try to solve code "katas" with it when it seems to fit well. But it has been my desktop calculator for everyday calculations for some time.

[robpike]

By the way, your id operator is very clear, but one conventional and slightly more obscure method is

3 3 rho 1 0 0 0

which is also

3 3 rho 4 take 1

and can be expressed several ways; here's one that doesn't need parentheses:

op id n = n n rho 1, n take 0

[fzipp]

By the way, your id operator is very clear

I learned it from this 1975 APL demonstration.

op id n = n n rho 1, n take 0

Thank you, this is clearly more efficient, as indicated by )cpu.

[robpike]

I have now pushed an implementation to GitHub. Commit 6e13ebc

[fzipp]

Thank you, that looks great! I'm closing the discussion, even though it's just a "Show and tell," as there is likely no further need for discussion on this matter.