[Numpy-discussion] Silent Broadcasting considered harmful

Sturla Molden sturla.molden at gmail.com
Mon Feb 9 03:47:26 EST 2015


On 09/02/15 08:34, Stefan Reiterer wrote:

> So maybe the better way would be not to add warnings to braodcasting
> operations, but to overhaul the matrix class
> to make it more attractive for numerical linear algebra(?)

I think you underestimate the amount of programming this would take. 
Take an arch LA trait of Matlab, the backslash operator.

If you have an expression like a \ b it is evaluated as solve(a,b). But 
how should you solve this? LU, QR, Cholesky, SVD? Exact? Least-squares? 
Surely that depends on the context and the array.

Matlab's backslash operators amounts to about 100,000 LOC.

Let's say we add attributes to the Matrix class to symbolically express 
things like inversion, symmetric matrix, triangular matrix, tridiagonal 
matrix, etc.

Say that you know a is symmetric and positive definite and b 
rectangular. How would you solve this most efficiently with BLAS and LAPACK?

    (a**-1) * b

In SciPy notation this is cho_solve(cho_factor(A),B). But what if you 
know that B has a certain structure?

And what if the arrays are in C order instead of Fortran order? Is there 
a way to avoid copy-and-transpose? NumPy's dot operator does this, but 
LA solvers should too. Put that on top of kernels for evalutationg any 
kind of (a**-1) * b expressions.

But what if the orders are reversed?

    b * (a**-1)

Those will also need special casings.

A LA package would need specialized code for any of these thinkable 
combinations, and pick the best one in a fly. That is what Matlab's 
backslash operator can do. NumPy cannot. But it is not difficult if you 
want to sit down and program it. If would just take you about 100,000 to 
half a million LOC. And chances are you effort will not even be appreciated.


Sturla






More information about the NumPy-Discussion mailing list