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

[Bob Dowling]

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.


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)