[Lang] Python-scope matrix/vector operations

[yuanming-hu]

Concisely describe the proposed feature
Currently, the Matrix class operations (e.g. matrix add/multiply/scale) can happen only in Taichi kernels. In some cases having these work in Python-scope is necessary. For example,
https://github.com/taichi-dev/taichi/pull/1006/files#diff-9ce279666131a93e7f2be366688b2fa1R52-R67

Without Python-scope matrix operation supports users have to fall back on numpy and the conversion between Taichi and numpy can be troublesome.

Describe the solution you'd like (if any)
Dispatch matrix operations according to the scope. Keep the old behavior if the operation happens in Taichi-scope, but evaluate the results immediately in Python-scope.

To decide the scope, use ti.get_runtime().inside_kernel. E.g.,

taichi/python/taichi/lang/expr.py

Line 143 in 6dce0c4

assert not ti.get_runtime().inside_kernel

[k-ye]

FWIW, when I created #1006 , I didn't intend to do the computation via ti.Matrix in the beginning. The reason is that the matrices in that example are static, so I feel like it's fine to use numpy to do the computation, then copy the data to ti.Matrix. My thought is that numpy's linalg ecosystem is much more mature, so if the matrix can be constructed once, we can just rely on numpy, and Taichi doesn't have to re-invent them... So actually I think it's fine to just support the common matrix ops in the Taichi scope.

users have to fall back on numpy and the conversion between Taichi and numpy can be troublesome.

I do agree with this. The current code looks messy for doing this conversion:

taichi/examples/cornell_box.py

Lines 65 to 67 in a1b397a

	m_inv = [list(r) for r in m_inv]
	m_inv_t = [list(r) for r in m_inv_t]
	return ti.Matrix(m_inv), ti.Matrix(m_inv_t)

So I feel like we can start by providing ti.Matrix.from_numpy() and ti.Matrix.to_numpy()? Excuse for my ignorance, i didn't know we already have these helpers...

[rexwangcc]

If we still want to do this after the discussion above, I'd like to take a look

[yuanming-hu]

My thought is that numpy's linalg ecosystem is much more mature, so if the matrix can be constructed once, we can just rely on numpy, and Taichi doesn't have to re-invent them... So actually I think it's fine to just support the common matrix ops in the Taichi scope.

I agree with this to some extend and I hate reinventing the wheels. However, a couple of users did complain about something like why can't we just do Taichi matrix inversion in Python scope? I feel like supporting these in native Taichi (instead of conversion to numpy) will help those people (especially for those who are not familiar with numpy).

Therefore I think it still makes sense to support Taichi matrix operations in Python scope, but we should start small: @rexwangcc if you are interested in this, how about just implementing matrix +, - and @ (matmul) as the first step, and let's see if there's any additional issue you encounter during development?

[rexwangcc]

Thanks! This sounds great to me, I personally also feel convenient to be able to do matrix operations in Python scope, and starting with +, -, @ should help me gradually understand how the frontend is implemented right now.

[rexwangcc]

Hi! I have some questions as I start to look into this.

It seems the instantiation of ti.Matrix is a bit messy (it also seems #1046 will simplify things a lot) and therefore caused some confusion of mine around the numpy io, since the first solution for @ came into my mind is simpy evaluate m1.to_numpy() @ m2.to_numpy(), I did some quick checks:

m1 = ti.Matrix(np.eye(n) * 2)
print(m1.to_numpy())

# File "taichi/python/taichi/lang/matrix.py", line 534, in to_numpy
#     ret = np.empty(self.loop_range().shape() + dim_ext,
# TypeError: 'tuple' object is not callable

m1 = ti.Matrix([[1, 2], [3, 4]])
print(m1.to_numpy())

# File "/home/rex/Code/GithubProjects/taichi-dev/taichi/python/taichi/lang/matrix.py", line 534, in to_numpy
#    ret = np.empty(self.loop_range().shape() + dim_ext,
# AttributeError: 'int' object has no attribute 'shape'

m1 = ti.Matrix(n, n, ti.f32, shape=(1, 1))
print(m1.to_numpy())
# worked as expected

It seems the underlying loop_range() which simply gets the first item of the entries is not working for the Matrix instance created from np array. (please correct me if I misunderstood anything here)

[yuanming-hu]

My feeling is that to_numpy() works for tensors of matrices but doesn't work for a single matrix. We will have to implement that in the future...

[rexwangcc]

Thanks, that makes sense to me. I guess at this point I may either implement that in order to use it for @ or simply relying on the values of self.entries, self.n, self.m to compute the matmul (which can be much slower than the numpy computation tho)

Actually let me go back and read the doc about external array again...