[Numpy-discussion] Confused by spec of numpy.linalg.solve

[Sebastian Berg]

On Di, 2014-04-01 at 15:31 +0100, Bob Dowling wrote:
> Versions:
> 
> >>> sys.version
> '3.3.2 (default, Mar  5 2014, 08:21:05) \n[GCC 4.8.2 20131212 (Red Hat
> 4.8.2-7)]'
> 
> >>> numpy.__version__
> '1.8.0'
> 
> 
> 
> Problem:
> 
> I'm trying to unpick the shape requirements of numpy.linalg.solve().
> The help text says:
> 
> solve(a, b) -
>      a : (..., M, M) array_like
>          Coefficient matrix.
>      b : {(..., M,), (..., M, K)}, array_like
>          Ordinate or "dependent variable" values.
> 
> It's the requirements on "b" that are giving me grief.  My read of the
> help text is that "b" must have a shape with either its final axis or
> its penultimate axis equal to M in size.  Which axis the matrix
> contraction is along depends on the size of the final axis of "b".
> 
> 
> So, according to my reading, if "b" has shape (6,3) then the first
> choice, "(..., M,)", is invoked but if "a" has shape (3,3) and "b" has
> shape (3,6) then the second choice, "(..., M, K)", is invoked.  I find
> this weird, but I've dealt with (much) weirder.
> 

I bet the documentation needs some more info there (if you have time,
please write a pull request). If you look at the code (that part is just
python code), you will see what really happens.

If `a` has exactly one dimension more then `b`, the first case is used.
Otherwise (..., M, K) is used instead. To make sure you always get the
expected result, it may be best to make sure that the number of
broadcasting (...) dimensions of `a` and `b` are identical (I am not
sure if you expect this to be the case or not). The shape itself does
not matter, only the (relative) number of dimensions does for the
decision which of the two signatures is used.

In other words, since you do not use `...` your examples always use the
(M, K) logic.

 - Sebastian

> 
> However, this is not what I see.  When "b" has shape (3,6) everything
> goes as expected.  When "b" has shape (6,3) I get an error message that
> 6 is not equal to 3:
> 
> > ValueError: solve: Operand 1 has a mismatch in its core dimension 0,
> > with gufunc signature (m,m),(m,n)->(m,n) (size 6 is different from 3)
> 
> 
> 
> Obviously my reading is incorrect.  Can somebody elucidate for me
> exactly what the requirements are on the shape of "b"?
> 
> 
> 
> Example code:
> 
> import numpy
> import numpy.linalg
> 
> # Works.
> M = numpy.array([
>      [1.0,     1.0/2.0, 1.0/3.0],
>      [1.0/2.0, 1.0/3.0, 1.0/4.0],
>      [1.0/3.0, 1.0/4.0, 1.0/5.0]
>      ] )
> 
> yy1 = numpy.array([
>      [1.0, 0.0, 0.0],
>      [0.0, 1.0, 0.0],
>      [0.0, 0.0, 1.0]
>      ])
> print(yy1.shape)
> xx1 = numpy.linalg.solve(M, yy1)
> print(xx1)
> 
> # Works too.
> yy2 = numpy.array([
>      [1.0, 0.0, 0.0, 1.0, 0.0, 0.0],
>      [0.0, 1.0, 0.0, 0.0, 1.0, 0.0],
>      [0.0, 0.0, 1.0, 0.0, 0.0, 1.0]
>      ])
> print(yy2.shape)
> xx2 = numpy.linalg.solve(M, yy2)
> print(xx2)
> 
> # Fails.
> yy3 = numpy.array([
>      [1.0, 0.0, 0.0],
>      [0.0, 1.0, 0.0],
>      [0.0, 0.0, 1.0],
>      [1.0, 0.0, 0.0],
>      [0.0, 1.0, 0.0],
>      [0.0, 0.0, 1.0]
>      ])
> print(yy3.shape)
> xx3 = numpy.linalg.solve(M, yy3)
> print(xx3)
> 
> 
> 
> 
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>