[Numpy-discussion] Fwd: Backslash operator A\b and np/sp.linalg.solve

[CJ Carey]

I agree that this seems more like a scipy feature than a numpy feature.

Users with structured matrices often use a sparse matrix format, though the
API for using them in solvers could use some work. (I have a
work-in-progress PR along those lines here:
https://github.com/scipy/scipy/pull/6331)

Perhaps this polyalgorithm approach could be used to dispatch sparse
matrices to the appropriate solver, while optionally checking dense
matrices for structure before dispatching them as well. Usage might look
like:

# if A is sparse, use scipy.sparse.linalg.solve, otherwise use
scipy.linalg.solve
scipy.linalg.generic_solve(A, b)

# converts A to banded representation and calls scipy.linalg.solveh_banded,
regardless if A is sparse or dense
scipy.linalg.generic_solve(A, b, symmetric=True, banded=(-5, 5))

# runs possibly-expensive checks, then dispatches to the appropriate solver
scipy.linalg.generic_solve(A, b, detect_structure=True)


(I'm not advocating for "generic_solve" as the final name, I just needed a
placeholder.)

On Tue, Jan 10, 2017 at 4:58 AM, Ilhan Polat <ilhanpolat at gmail.com> wrote:

> I've done some benchmarking and it seems that the packed storage comes
> with a runtime penalty which agrees with a few links I've found online
>
> https://blog.debroglie.net/2013/09/01/lapack-and-packed-storage/
> http://stackoverflow.com/questions/8941678/lapack-are-
> operations-on-packed-storage-matrices-faster
>
> The access of individual elements in packed stored matrices is expected to
> be more costly than in full storage, because of the more complicated
> indexing necessary. Hence, I am not sure if this justifies the absence just
> by having a dedicated solver for a prescribed structure.
>
> Existence of these polyalgorithms in matlab and not having in lapack
> should not imply FORTRAN users always know the structure in their matrices.
> I will also ask in LAPACK message board about this for some context.
>
> But thanks tough. As usual there is more to it than meets the eye
> probably,
> ilhan
>
>
>
>
> On Tue, Jan 10, 2017 at 4:16 AM, Robert Kern <robert.kern at gmail.com>
> wrote:
>
>> On Mon, Jan 9, 2017 at 7:10 PM, Ilhan Polat <ilhanpolat at gmail.com> wrote:
>> >
>> > Yes, that's precisely the case but when we know the structure we can
>> just choose the appropriate solver anyhow with a little bit of overhead.
>> What I mean is that, to my knowledge, FORTRAN routines for checking for
>> triangularness etc. are absent.
>>
>> I'm responding to that. The reason that they don't have those FORTRAN
>> routines for testing for structure inside of a generic dense matrix is that
>> in FORTRAN it's more natural (and efficient) to just use the explicit
>> packed structure and associated routines instead. You would only use a
>> generic dense matrix if you know that there isn't structure in the matrix.
>> So there are no routines for detecting that structure in generic dense
>> matrices.
>>
>> --
>> Robert Kern
>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> https://mail.scipy.org/mailman/listinfo/numpy-discussion
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20170110/cb7d8ade/attachment.html>