New result shaping logic fails for non-square A

[Technologicat]

In the updated spqr.qr_solve, the result is formatted by

    ## Return x with the same shape as b.
    numpy_x = cholmoddense2numpy( chol_x ).reshape( b.shape )

but this fails if A is not square. Consider A.shape = (m, n) (m equations, n unknowns), b.shape = (m, k) (m equations, k independent RHS).

Then, if k > 1, and we use the same "independent problems as columns" convention as in b, we should have x.shape = (n, k).

In the k = 1 case, I think that to keep with the spirit of NumPy, it makes the most sense to return a rank-1 array.

If this is the preferred way to fix this, the above could be replaced by e.g. (tested, works):

    ## Return x with the appropriate shape.
    ## The result is an n-by-k matrix if k > 1; for a single RHS, just a vector.
    B = numpy.atleast_2d(b)
    if B.shape[0] == 1  and  B.shape[1] > 1:
        B = B.T  # prefer column vector
#    m = A.shape[0] # == B.shape[0]
    n = A.shape[1]
    k = B.shape[1]
    shp = (n,k) if k > 1 else (n,)
    numpy_x = cholmoddense2numpy( chol_x ).reshape( shp )

but this solution is a bit verbose, and in the k > 1 case reshapes an array that is already in the correct shape.

Hence, I'm raising the issue for discussion: what is the best solution here?

And do we need to update also qr_solve_sparse to do something similar? (I think "no" since, IIRC, SciPy sparse matrices are always 2D.)

[yig]

Good catch. I think the dimension of the result should match the dimension of the input. So a rank-1 array in produces a rank-1 array out, but a rank-2 array in produces a rank-2 array out, even if the second dimension is 1. Here is my solution (tested, pushed):

I replaced:

numpy_x = cholmoddense2numpy( chol_x ).reshape( b.shape )

with:

x_shape = list( b.shape )
x_shape[0] = A.shape[1]
numpy_x = cholmoddense2numpy( chol_x ).reshape( x_shape )

[Technologicat]

Nice, compact solution. I see you also updated the test.

Preserving input rank is a good point. This is probably the best choice for making the behavior predictable, and very much in line with the spirit of NumPy.

I pulled the changes to my fork, too. Thanks!