Matrix reshaping

[ghost]

Depends on #3 if you think in column major terms (common vector algebra habit).

From Octave, MATLAB, FORTRAN:

>> help reshape
'reshape' is a built-in function from the file libinterp/corefcn/data.cc

 -- Built-in Function: reshape (A, M, N, ...)
 -- Built-in Function: reshape (A, [M N ...])
 -- Built-in Function: reshape (A, ..., [], ...)
 -- Built-in Function: reshape (A, SIZE)
     Return a matrix with the specified dimensions (M, N, ...) whose
     elements are taken from the matrix A.

     The elements of the matrix are accessed in column-major order (like
     Fortran arrays are stored).

     The following code demonstrates reshaping a 1x4 row vector into a
     2x2 square matrix.

          reshape ([1, 2, 3, 4], 2, 2)
                =>  1  3
                    2  4

     Note that the total number of elements in the original matrix
     ('prod (size (A))') must match the total number of elements in the
     new matrix ('prod ([M N ...])').

     A single dimension of the return matrix may be left unspecified and
     Octave will determine its size automatically.  An empty matrix ([])
     is used to flag the unspecified dimension.

     See also: resize, vec, postpad, cat, squeeze.

[tafia]

I think a row major reshaped MatrixSlice is a no-op. I am not familiar with the reshape function itself, can't the algorithms using it be adapted to a row-major version? (we must discourage users from using a slow function if possible ...)

[ghost]

Well, the zero-effort version is matrix.transpose().reshape() to get an equivalent, but as I understand it, in-place transposition is roughly O(MN) or worse, but I guess there are papers about an O(n) algorithm. I haven't read any of them, but "linear in the number of elements" takes you back to O(MN).

I don't know what rulinalg is doing right now, but I imagine that knowing the user's intention to do columnwise things would probably improve performance in that case.

[Andlon]

I think it makes sense to provide some support for reshaping. As @tafia points out, it's a no-op if the data layout used is the same (i.e. row-major), and this would be relatively easy to implement as far as I can tell.

However, let's just step back for a second and consider where reshaping is used. In my experience, reshape in MATLAB/numpy is used primarily because every operation in these languages must be done with whole data sets (i.e. "vectorized" in the sense of batch operations on a number of elements with each operation) in order to show any kind of efficiency, since elementwise access is very, very slow in these languages. For example, you might construct a matrix by cleverly constructing a few other kinds of matrices, combining them in some way, and reshaping the result to fit your needs.

In a language like Rust, where elementwise access is generally (barring SIMD) as fast as vector operations (in the MATLAB sense), and often faster due to the avoidance of temporary allocations of relatively big matrices, you can often construct your matrix directly in terms of its elements with no performance penalty. This greatly alleviates the need for this kind of functionality.

The above is based on my own experience and usage patterns. I'm interested to hear more about how reshape is used by others. In particular, I'd like to see examples where reshaping is considerably more convenient than direct construction!

An alternative to providing reshape functionality could be to provide documentation for how to rewrite reshape-based matrix construction in MATLAB/numpy to direct elementwise construction in rulinalg. I don't know (yet) which is more appropriate!

[AtheMathmo]

As @tafia said we should be able to provide an equivalent no-op reshape function. It would simply work in row-major instead of column-major order.

@bright-star I did play around with in-place transposition. The short of it is that even with a fairly smart algorithm you'd still expect it to be slower than an out-of-place transposition. Generally a strided transpose would be better (which we don't support) and in-place becomes especially valuable when we need to move the data around and it is very large.

And finally I do agree with @Andlon that I'm not sure of the need for reshape. This is probably more from my own experience and I'd also be happy to be convinced otherwise.

Regardless, I think this could make a valuable addition to the BaseMatrix trait. Something like:

fn reshape(self, rows: usize, cols: usize) -> Self { ... }

I'm not sure if we can actually reshape a MatrixSlice? If not we can just stick with Matrix.

[ghost]

Here's an immediate use case (scroll to the end to see the matrix equations, which are just Ax = b): https://bright-star.gitlab.io/electro/finite-difference/

What happens in that application is that you have a grid of points that directly describe a rectangular physical geometry (say, m rows and n columns). However, the equations you need form coefficients with a Toeplitz matrix in the order of the number of points, so, it's mn x mn. Since the coefficients are Toeplitz-symmetric, majority doesn't matter. However, you need to form the b vector from the known values in the physical geometry, by taking each of the columns and "pouring out" the matrix of knowns into a single mn x 1 vector.

A single call to reshape accomplishes this: bounds = reshape(mat,node_count,1);

Currently I'm using: let boundary_vec = Vector::<f64>::from(mat.transpose().into_vec());

---

There are more dramatic examples in tensor algebra, such as fiber products.

[AtheMathmo]

Ah I see - though I think we'd need a fair amount of work for this kind of reshape to be accessible. I think the code you have there is probably about as efficient as you can get without implementing some kind of column-major strided view (see #26 ). If we were to implement reshape I think that it would need to keep the matrix row-major for efficiency. I've avoided column-major alternatives so far as implementing them efficiently is a lot of work, sadly.

Of course if you transposed your algorithm so that you can "pour out" rows. You'd have a no-op here:

let boundary_vec = Vector::<f64>::from(mat.into_vec());

(Though this could be on a similar scale of complexity as supporting column-major in the first place!)

[ghost]

Right. You end up pushing the bubble of effort around, basically.

My personal opinion is that the application code should at least resemble the mathematics semantically. A big reason why people use Octave or MATLAB despite the proprietary cost and technical issues is the intuitiveness granted by commitment to vectors and matrices as first-class objects. The truth is that it's not great and the syntax is very error-prone, but there should also be a usability reward for a mathematical programmer to switch away from Octave or MATLAB. The safety and performance of Rust is great but when an implementation pipeline is "Write it in MATLAB, see if it works, then write it in C", it's not enough to convince.

R is a good example. It is also a vector language with great linalg bindings and precompiled libs and so on, but it looks like Javascript and has a bunch of functions for taking lambdas over different vector slices. So it has good usability and it's better than Octave or MATLAB for serious data analysis.

In light of all this, I think that there's no need to try accommodate column-major as an actual way of structuring the data, but if it were possible to write API functions that allowed the user to act as though the matrix were column-major without having to do extra work, everyone might be happy.

[Andlon]

Thanks for the example @bright-star! It's been a while since I worked with finite differences.

I essentially agree with your assessment that reshaping is a useful addition. I still believe it is rarely the best tool for the job, but I guess it's somewhat intuitive to use in some cases.

I also think that we should not worry too much about the efficiency of a possible implementation. If you reshape using a row-major order, great, it's a no-op. If you want to reshape using column-major order, you have to pay the price. However, compared to, for example, solving a linear system with the same matrix, this cost should be negligible.

I also note that in your example @bright-star, you want to reshape a matrix as a vector. This wouldn't actually fit with the proposed addition, since matrix.reshape() would return a matrix. I guess we'd need a separate method on the vector, i.e. vec.flatten() or so.

EDIT: vec.flatten() should probably rather have been mat.flatten(), in that case!

[ghost]

@Andlon That's absolutely correct. From the typed language perspective (and from a math perspective), reshaping into a vector and not explicitly handling it as one is a sloppy habit. I actually think that the Rust/rulinalg idiom here is perfect as-is.

As proof, here's what I actually did to form the constant vector for that FD problem so that .solve() would accept it:

let boundary_vec = Vector::<f64>::from(
    pot_mat.transpose()
           .into_vec()
           .iter()
           .map(|el| { el * -1.0})
           .collect::<Vec<f64>>());

[ghost]

I should also point out for the record (I'm sure you already know) that the column-major convention leads to left matrix multiplication and division for these associated problems, but if I remember correctly, row-major convention just changes all left-side operations to right-side. So correctness is not hard to verify. (Do I have that right?)

[Andlon]

I should also point out for the record (I'm sure you already know) that the column-major convention leads to left matrix multiplication and division for these associated problems, but if I remember correctly, row-major convention just changes all left-side operations to right-side. So correctness is not hard to verify. (Do I have that right?)

@bright-star: sorry, I'm not following. What exactly do you mean by associated problems?

[ghost]

Sorry, I meant linear algebra problems. If a matrix is considered a collection of row vectors instead of column vectors, then b = xA, where x and b are row vectors, is the same equation as its column vector equivalent transpose(A)transpose(x) = transpose(b). And this can be applied mechanically everywhere. Isn't that the case?

[tafia]

Solving a linear equation or its transposed version is the same cost (+/- a no-op reshaping to convert the "row vector" to a "column vector").

But I understand that people not familiar with these technical details should be able to get the job done.

We could probably have 2 reshape functions and state in the docs that the row major version is a no-op.

On 25 Oct 2016 04:42, "Andreas Borgen Longva" notifications@github.com wrote:

Thanks for the example @bright-star https://github.com/bright-star! It's been a while since I worked with finite differences.

I essentially agree with your assessment that reshaping is a useful addition. I still believe it is rarely the best tool for the job, but I guess it's somewhat intuitive to use in some cases.

I also think that we should not worry too much about the efficiency of a possible implementation. If you reshape using a row-major order, great, it's a no-op. If you want to reshape using column-major order, you have to pay the price. However, compared to, for example, solving a linear system with the same matrix, this cost should be negligible.

I also note that in your example @bright-star https://github.com/bright-star, you want to reshape a matrix as a vector. This wouldn't actually fit with the proposed addition, since matrix.reshape() would return a matrix. I guess we'd need a separate method on the vector, i.e. vec.flatten() or so.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/AtheMathmo/rulinalg/issues/72#issuecomment-255860063, or mute the thread https://github.com/notifications/unsubscribe-auth/AHAszjyyZk3dcPYyrfHJD1cV-_D474tiks5q3RgxgaJpZM4KeU0R .