[Numpy-discussion] is_triangular, is_diagonal, is_symmetric et al. in NumPy or SciPy linalg

Fri Jul 2 07:19:24 EDT 2021

Ah right. So two things, the original reason f9r this question is because I
can't decide in https://github.com/scipy/scipy/pull/12824 whether others
would also benefit from quick structure determination.

I can keep it private function or we can put them some misc or lib folder
so all can use. Say there is a special method for triangular matrices but
you can't guarantee the structure so you can quickly check for it. At worst
O(n**2) complexity for diagonal arrays and almost O(2n) for full arrays
makes it quite appealing.

But then again maybe NumPy is a better place since probably it will be
faster to have this in pure C with the right headers and without the extra
Cython overhead.

Funny you mention the container idea. This is precisely what I'm doing in
PR mentioned above (I'll push when I'm done). I stole the idea from Tim
Davis himself in a Julia discussion for keeping the factorization as an
attribute to be used later if need be. So yes it makes a lot of sense
Sparse or not.

On Wed, 30 Jun 2021, 19:14 Evgeni Burovski, <evgeny.burovskiy at gmail.com>
wrote:

> Hi Ilhan,
>
> Overall I think something like this would be great. However, I wonder
> if you considered having a specialized container with a structure tag
> instead of trying to discover the structure. If it's a container, it
> can neatly wrap various lapack storage schemes and dispatch to an
> appropriate lapack functionality. Possibly even sparse storage
> schemes. And it seems a bit more robust than trying to discover the
> structure (e.g. what about off-band elements of  \sim 1e-16 etc).
>
> The next question is of course if this should live in scipy/numpy
> .linalg or as a separate repo, at least for some time (maybe in the
> scipy organization?). So that it can iterate faster, among other
> things.
> (I'd be interested in contributing FWIW)
>
> Cheers,
>
> Evgeni
>
>
> On Wed, Jun 30, 2021 at 1:22 AM Ilhan Polat <ilhanpolat at gmail.com> wrote:
> >
> > Dear all,
> >
> > I'm writing some helper Cythpm functions for scipy.linalg which is kinda
> performant and usable. And there is still quite some wiggle room for more.
> >
> > In many linalg routines there is a lot of performance benefit if the
> structure can be discovered in a cheap and reliable way at the outset. For
> example if symmetric then eig can delegate to eigh or if triangular then
> triangular solvers can be used in linalg.solve and lstsq so forth
> >
> > Here is the Cythonized version for Jupyter notebook to paste to discover
> the lower/upper bandwidth of square array A that competes well with A != 0
> just to use some low level function (note the latter returns an array hence
> more cost is involved) There is a higher level supervisor function that
> checks C-contiguousness otherwise specializes to different versions of it
> >
> > Initial cell
> >
> > %load_ext Cython
> > %load_ext line_profiler
> > import cython
> > import line_profiler
> >
> > Then another cell
> >
> > %%cython
> > # cython: language_level=3
> > # cython: linetrace=True
> > # cython: binding = True
> > # distutils: define_macros=CYTHON_TRACE=1
> > # distutils: define_macros=CYTHON_TRACE_NOGIL=1
> >
> > cimport cython
> > cimport numpy as cnp
> > import numpy as np
> > import line_profiler
> > ctypedef fused np_numeric_t:
> >     cnp.int8_t
> >     cnp.int16_t
> >     cnp.int32_t
> >     cnp.int64_t
> >     cnp.uint8_t
> >     cnp.uint16_t
> >     cnp.uint32_t
> >     cnp.uint64_t
> >     cnp.float32_t
> >     cnp.float64_t
> >     cnp.complex64_t
> >     cnp.complex128_t
> >     cnp.int_t
> >     cnp.long_t
> >     cnp.longlong_t
> >     cnp.uint_t
> >     cnp.ulong_t
> >     cnp.ulonglong_t
> >     cnp.intp_t
> >     cnp.uintp_t
> >     cnp.float_t
> >     cnp.double_t
> >     cnp.longdouble_t
> >
> >
> > @cython.linetrace(True)
> > @cython.initializedcheck(False)
> > @cython.boundscheck(False)
> > @cython.wraparound(False)
> > cpdef inline (int, int) band_check_internal(np_numeric_t[:, ::1]A):
> >     cdef Py_ssize_t n = A.shape[0], lower_band = 0, upper_band = 0, r, c
> >     cdef np_numeric_t zero = 0
> >
> >     for r in xrange(n):
> >         # Only bother if outside the existing band:
> >         for c in xrange(r-lower_band):
> >             if A[r, c] != zero:
> >                 lower_band = r - c
> >                 break
> >
> >         for c in xrange(n - 1, r + upper_band, -1):
> >             if A[r, c] != zero:
> >                 upper_band = c - r
> >                 break
> >
> >     return lower_band, upper_band
> >
> > Final cell for use-case ---------------
> >
> > # Make arbitrary lower-banded array
> > n = 50 # array size
> > k = 3 # k'th subdiagonal
> > R = np.zeros([n, n], dtype=np.float32)
> > R[[x for x in range(n)], [x for x in range(n)]] = 1
> > R[[x for x in range(n-1)], [x for x in range(1,n)]] = 1
> > R[[x for x in range(1,n)], [x for x in range(n-1)]] = 1
> > R[[x for x in range(k,n)], [x for x in range(n-k)]] = 2
> >
> > Some very haphazardly put together metrics
> >
> > %timeit band_check_internal(R)
> > 2.59 µs ± 84.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops
> each)
> >
> > %timeit np.linalg.solve(R, zzz)
> > 824 µs ± 6.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
> >
> > %timeit R != 0.
> > 1.65 µs ± 43.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops
> each)
> >
> > So the worst case cost is negligible in general (note that the given
> code is slower as it uses the fused type however if I go with tempita
> standalone version is faster)
> >
> > Two questions:
> >
> > 1) This is missing np.half/float16 functionality since any arithmetic
> with float16 is might not be reliable including nonzero check. IS it safe
> to view it as np.uint16 and use that specialization? I'm not sure about the
> sign bit hence the question. I can leave this out since almost all linalg
> suite rejects this datatype due to well-known lack of supprt.
> >
> > 2) Should this be in NumPy or SciPy linalg? It is quite relevant to be
> on SciPy but then again this stuff is purely about array structures. But if
> the opinion is for NumPy then I would need a volunteer because NumPy
> codebase flies way above my head.
> >
> >
> > All feedback welcome
> >
> > Best
> > ilhan
> >
> > _______________________________________________
> > NumPy-Discussion mailing list
> > NumPy-Discussion at python.org
> > https://mail.python.org/mailman/listinfo/numpy-discussion
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at python.org
> https://mail.python.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/numpy-discussion/attachments/20210702/aa3d9a21/attachment-0001.html>