Linear Algebra routines as generalized ufuncs #2954

Closed
wants to merge 77 commits into
from

Conversation

Projects
None yet
Contributor

ovillellas commented Jan 31, 2013

This pull request contains a module with some lineal algebra functions implemented as generalized ufuncs. It is implemented using LAPACK/BLAS, and it is configured to compile using the same LAPACK/BLAS configuration as numpy.linalg. It also includes some extras, like some fused arithmetic ufuncs (multiply_add and similar ones).

More information can be found in gufuncs_linalg_contents.rst

mwiebe and others added some commits Oct 19, 2012

@mwiebe mwiebe BUG: Fix generalized ufunc so repeating labels in one input/output is ok 1613e09
@mwiebe mwiebe BUG: Missed case of broadcasting onto core dimensions in gufunc be911e8
@mwiebe mwiebe BUG: Make sure the nditer doesn't complain when 0 broadcast dims 9a70cfd
@mwiebe mwiebe TST: Add test for gufunc scalar case b5c7733
@mwiebe mwiebe BUG: Fix bug in gufunc scalar case 12be5b7
@mwiebe mwiebe TST: Test for a generalized ufunc bug, for zero-sized inputs fa9dbef
@mwiebe mwiebe BUG: Fix for generalized ufunc zero-sized input case cac3de5
@ovillellas ovillellas created a new module to hold linalg ufuncs. 153e02e
@ovillellas ovillellas inner1d and mat_mult implemented using blas. 06613da
@ovillellas ovillellas refactored some code, make it cleaner overall and ready to reuse some…
… code from the matrix_multiply in other gufuncs
0ab4a14
@ovillellas ovillellas det and slogdet working e6d08ea
@ovillellas ovillellas eigh and eigvalsh working a3baf1c
@ovillellas ovillellas solve gufunc working e53b3d3
@ovillellas ovillellas solve1 and inv working 6ad2738
@ovillellas ovillellas fixed possible issues with BLAS _copy (0 is not a valid incx/incy val…
…ue and the functions are not guaranteed to work with them).

also got cholesky working.
1276782
@ovillellas ovillellas working eig and eigvals priority 2 functions. ee4138d
@ovillellas ovillellas svd implemented. Single output working. Multiple options not function…
…al due to a bug in the harness.
b1c3e09
@ovillellas ovillellas modified the code so it just used external definitions of blas/lapack…
… functions (as is made in linalg module). Changed some code so that it uses BLAS instead of cblas (the fortran functions) also in line with what it is done on linalg. Modified the matrix multiply code, made it simpler and adapted to use blas (it was using cblas with CblasRowMajor, that is not available in the fortran function.
c229ba0
@ovillellas ovillellas lapack_lite for builds of umath_linalg without an optimized lapack in…
… the system.
4b5f777
@ovillellas ovillellas added some single precision functions to f2c_lite.c that were missing…
… and needed by out library. It seems to work now on Linux.
b29fa3c
@ovillellas ovillellas added plenty of simple functions (quadratic_form plus all the "inspir…
…ed from PDL" ufuncs). Only missing from "inspired from PDL" is matrix_multiply3.
56d70ff
@ovillellas ovillellas added information about the contents of umath-linalg module 3cde140
@ovillellas ovillellas fixed gufuncs so that they use the proper signature (mwiebe fix present) f26444d
@ovillellas ovillellas fixed a warning in f2c_lite.c for umath/lapack_lite
added chosolve and poinv, they need testing.
839e597
@ovillellas ovillellas poinv and chosolve working. Rebuilt lapack_lite to support them. Used…
… also a patched f2c to remove warnings.
4f2a508
@ovillellas ovillellas updated umath_linalg_content.txt b7d0f7b
@ovillellas ovillellas wrote a wrapper module for umath_linalg. Named gufuncs_linalg (in pyt…
…hon).
6100e8a
@ovillellas ovillellas first iteration with tests. Incomplete and some failing. Just a start…
…. Some bugs already fixed.
fdeafbd
@ovillellas ovillellas work on tests and related fixes. Getting things in shape to commit to…
… de-shaw patch
3bfb914
@ovillellas ovillellas removed some wrappers that weren't needed with the harness fix, just …
…changed to assignment
c4a2cff
@ovillellas ovillellas modified umath_linalg_content.txt to reflect changes. 23a6380
@ovillellas ovillellas updated the umath_linalg_content.txt adding a mention to the wrapper …
…code.
0735c95
@ovillellas ovillellas fixed bug in matrix_multiply when using cdoubles 6e0ea34
@ovillellas ovillellas fixed the problem in eigvals (apparently) 250786b
@ovillellas ovillellas work in progress: proper tests for eig. 8565a63
@ovillellas ovillellas added tests for ufuncs in gufuncs_linalg (the ones based on pdl). Add…
…ed multiply4 in the wrapper, as it was missing
f9372ab
@ovillellas ovillellas reverted matrix_dot in umath_gufuncs to matrix_multiply.
added some type-tests on test_gufuncs_linalg.py
d8ed21d
@ovillellas ovillellas updated documentation a5e5833
@ovillellas ovillellas compile fix on linux/mac f2c81b4
@ovillellas ovillellas Merge branch 'gufunc-fix' into gufunc-linalg-stable c25e0b7
@ovillellas ovillellas BLD: Windows build fixes + some tabs removed a4f540d
@ovillellas ovillellas STY: made sure that split strings had \ at the end 6a09b5c
@ovillellas ovillellas updated api version, as one merge changed it. e7b6a2b
@ovillellas ovillellas fixed testdet test. It failed due to eigvails failing in single preci…
…sion and notifying the failure as nans.
05d0146
@ovillellas ovillellas TST: fixed test for gufuncs_linalg Det 32b0617
@ovillellas ovillellas ENH: cholesky handling of _potrf failures (set result to nan) 22b163c
@ovillellas ovillellas ENH: eigh, eigvalsh set result to nan on LAPACK error (_ssyevd, _heevd) 6e1b48c
@ovillellas ovillellas ENH: solve sets result to nan on LAPACK error (_gesv) 69f6464
@ovillellas ovillellas ENH: inv sets result to nan on LAPACK error (_gesv) 48c4ab4
@ovillellas ovillellas ENH: svd sets results to nan on LAPACK error (_gesdd) a5a1b83
@ovillellas ovillellas ENH: chosolve sets result to nan on LAPACK error (_potrf, _potrs) 40d3746
@ovillellas ovillellas DOC: Added docstring for eigh 01186ca
@ovillellas ovillellas DOC: Added docstring to eigvalsh d7b2a58
@ovillellas ovillellas DOC: added docstring to solve a6c51c6
@ovillellas ovillellas DOC: added docstring for svd 89083d4
@ovillellas ovillellas DOC: Added docstring to chosolve e538255
@ovillellas ovillellas DOC: added docstring for poinv f008fdb
@ovillellas ovillellas DOC: Added notes on error handling. fa1c5f8
@ovillellas ovillellas MAINT: renamed umath_linalg module to _umath_linalg as it is internal. 3a970dc
@ovillellas ovillellas MAINT: renamed the file describing the gufuncs_linalg module 6cf6586
@ovillellas ovillellas MAINT: Rewrote the gufuncs_linalg_contents as a rst file and updated it. b686a76

@rkern rkern commented on an outdated diff Jan 31, 2013

numpy/core/src/umath/gufuncs_linalg.py
+"""Linear Algebra functions implemented as gufuncs, so they can be broadcast.
+
+Notes
+-----
+This module contains functionality that could be found in the linalg module,
+but in a gufunc-like way. This allows the use of vectorization and broadcasting
+on the operands.
+
+This module itself provides wrappers to kernels written as numpy
+generalized-ufuncs that perform the heavy-lifting of computation. The wrappers
+exist in order to provide a sane interface, like providing keyword arguments in
+line with the ones used by linalg, or just to automatically select the
+appropriate kernel depending on the parameters. All wrappers forward the keyword
+parameters to the underlying generalized ufunc (the kernel).
+
+The functions are intended to be used on arrays of functions. For those
@rkern

rkern Jan 31, 2013

Member

Did you mean "arrays of matrices" instead of "arrays of functions"?

@rkern rkern and 1 other commented on an outdated diff Jan 31, 2013

numpy/core/src/umath/gufuncs_linalg.py
+
+ Implemented for types single, double, csingle and cdouble. Numpy
+ conversion rules apply.
+
+ See Also
+ --------
+ multiply3 : element-wise three-way multiplication.
+ multiply3_add : element-wise three-way multiplication and addition.
+ multiply_add : element-wise multiply-add.
+ multiply_add2 : element-wise multiplication with two additions.
+ multiply4 : element-wise four-way multiplication
+ multiply4_add : element-wise four-way multiplication and addition,
+
+ Examples
+ --------
+ <Some example in doctest format>
@rkern

rkern Jan 31, 2013

Member

I would recommend either filling this out or removing the section before merging the PR.

@ovillellas

ovillellas Jan 31, 2013

Contributor

On Thu, Jan 31, 2013 at 5:52 PM, Robert Kern notifications@github.comwrote:

In numpy/core/src/umath/gufuncs_linalg.py:

  • Implemented for types single, double, csingle and cdouble. Numpy
  • conversion rules apply.
  • See Also

  • multiply3 : element-wise three-way multiplication.
  • multiply3_add : element-wise three-way multiplication and addition.
  • multiply_add : element-wise multiply-add.
  • multiply_add2 : element-wise multiplication with two additions.
  • multiply4 : element-wise four-way multiplication
  • multiply4_add : element-wise four-way multiplication and addition,
  • Examples

I would recommend either filling this out or removing the section before
merging the PR.

In fact the idea is to fill it before the merge. I also imagine that there
will be feedback anyways, so I wanted to start the pull request as I
imagine there will be changes to be made overall.

That's also why I split the two pull-request. I imagine that a new module
will require more changes than a fix :)


Reply to this email directly or view it on GitHubhttps://github.com/numpy/numpy/pull/2954/files#r2846449.

@rkern

rkern Jan 31, 2013

Member

Great. I just wanted to make sure it was explicitly on the TODO list. :-)

@njsmith njsmith commented on the diff Jan 31, 2013

numpy/core/setup.py
@@ -919,6 +919,38 @@ def get_dotblas_sources(ext, build_dir):
sources = [join('src','umath', 'umath_tests.c.src')])
#######################################################################
+ # umath_linalg module #
+ #######################################################################
+
+ lapack_info = get_info('lapack_opt', 0)
+ def get_lapack_lite_sources(ext, build_dir):
+ if not lapack_info:
+ print("### Warning: Using unoptimized lapack ###")
+ return ext.depends[:-1]
+ else:
+ if sys.platform=='win32':
+ print("### Warning: python_xerbla.c is disabled ###")
@njsmith

njsmith Jan 31, 2013

Owner
  • What is python_xerbla.c?
  • What effect does disabling it have?
  • Why is it disabled?
@ovillellas

ovillellas Feb 1, 2013

Contributor

In LAPACK xerbla is an error handling function. LAPACK routines print out their errors using that function. python_xerbla is a version of the function friendly with Python.

I had to make a version of lapack_lite, following the instructions present in lapack_lite inside the linalg package. I needed to do so as I use some routines not present on the current lapack_lite (notably the single precision ones and some extra).

The build configuration is inspired by the build configuration of linalg (i.e. ruthlessly copied :) ).

@njsmith njsmith and 1 other commented on an outdated diff Jan 31, 2013

numpy/core/src/umath/gufuncs_linalg.py
+-----
+This module contains functionality that could be found in the linalg module,
+but in a gufunc-like way. This allows the use of vectorization and broadcasting
+on the operands.
+
+This module itself provides wrappers to kernels written as numpy
+generalized-ufuncs that perform the heavy-lifting of computation. The wrappers
+exist in order to provide a sane interface, like providing keyword arguments in
+line with the ones used by linalg, or just to automatically select the
+appropriate kernel depending on the parameters. All wrappers forward the keyword
+parameters to the underlying generalized ufunc (the kernel).
+
+The functions are intended to be used on arrays of functions. For those
+functions where a result may not be possible to obtain (like the inverse of
+a matrix that is not invertible) no exception is raised, but the results for
+the elements involved are set to NaN.
@njsmith

njsmith Jan 31, 2013

Owner

Surely this should be determined by the value of np.seterr? Why have special rules here?

@ovillellas

ovillellas Feb 1, 2013

Contributor

May be I should reword this. The cases were the results are set to NaNs are the cases where linalg would raise a LinAlgError. That is, a method not converging, so no valid result would be given. np.seterr deals with floating point unit errors AFAIK (overflow, underflow and the like). Maybe make sure that when the code that fills with nans gets called that underlying error handling gets a nan (with feraiseexcept(FE_INVALID) maybe? )

The idea of this submodule is being able to call the kernels on an array of problems. Raise an exception when there was an error in one element did not seem appropriate.

@njsmith njsmith commented on an outdated diff Jan 31, 2013

numpy/core/src/umath/gufuncs_linalg.py
+
+ Multiplying on the left by an orthogonal matrix, `Q`, and on the
+ right by `Q.T` (the transpose of `Q` preserves the eigenvalues of
+ the original matrix
+
+ >>> x = np.random.random()
+ >>> Q = np.array([[np.cos(x), -np.sin(x)], [np.sin(x), np.cos(x)]])
+ >>> A = matrix_multiply(Q, D)
+ >>> A = matrix_multiply(A, Q.T)
+ >>> eigvals(A)
+ array([-1.+0.j, 1.+0.j])
+
+ """
+ return _impl.eigvals(a, **kwargs)
+
+def quadratic_form(u,Q,v, **kwargs):

@njsmith njsmith commented on the diff Jan 31, 2013

numpy/core/src/umath/gufuncs_linalg_contents.rst
+- An uniform parameter is needed to configure the way the computation is
+performed (like UPLO in the functions working on symmetric/hermitian matrices)
+- svd, where it was impossible to map the function to a gufunc signature.
+
+In the case of uniform parameters like UPLO, there are two separate entry points
+in the C module that imply either 'U' or 'L'. The wrapper just selects the
+kernel to use by checking the appropriate keyword parameter. This way a
+function interface similar to numpy.linalg can be kept.
+
+In the case of SVD not only there were problems with the support of keyword
+arguments. There was the added problem of the signature system not being able
+to cope with the needs of this functions. Just for the singular values a
+a signature like (m,n)->(min(m,n)) was needed. This has been worked around by
+implementing different kernels for the cases where min(m,n) == m and where
+min(m,n) == n. The wrapper code automatically selects the appropriate one.
+
@njsmith

njsmith Jan 31, 2013

Owner

Would it be possible to extend the core gufunc system to handle some or all of these problems? I always worry when we start making workarounds like this that we're just digging ourselves deeper into a problem we'll have to face sooner or later. In particular wanting to pass extra arguments to gufuncs in general seems pretty reasonable... or can we support a signature like (m,n)->(k) where the gufunc gets handed m and n and gets to compute k, before ever calling the actual loop?

@ovillellas

ovillellas Feb 1, 2013

Contributor

In fact we discussed that internally. It is beyond the scope of the module to change that, though.
It would certainly be interesting to give the gufunc implementer the ability to complete (and maybe also verify) the signature. That would allow supporting (m,n)->(k) in a way that computes k from m and n (the case of min(m,n) could be handled). It could also allow having that k as 2*n, or even things like checking that m <= n. IMO, that would be great (and could be implemented as optional).
It would also be nice to have a setup and a teardown function. This would allow allocating buffers and deallocating them once for all the execution. Right now, when running several dimensions that cannot be set as a single lineal iteration in the kernel that setup/teardown needs to be repeated in each call to the gufunc function. In this module per-element buffers are needed for intermediate results, that get reused. It is not much of a problem as it is just few allocations and they are not repeated per element, though.

@jaimefrio

jaimefrio Feb 22, 2014

Member

It's not perfect, but if you define a signature like (m,n)->k and call the kernel passing in an out array, then the value of k passed to the kernel will come from the shape of the out array. So you could have the shape calculation logic in the wrapper, create the output array there, and always call the kernel with an out parameter.

@njsmith njsmith commented on the diff Jan 31, 2013

numpy/core/src/umath/ufunc_object.c
@@ -71,6 +71,12 @@
static int
_does_loop_use_arrays(void *data);
+static int
+assign_reduce_identity_zero(PyArrayObject *result);
+
+static int
+assign_reduce_identity_one(PyArrayObject *result);
+
@njsmith

njsmith Jan 31, 2013

Owner

I guess this code will disappear from this PR's diff after we merge the other one...

@ovillellas

ovillellas Feb 4, 2013

Contributor

Yes, internally the branch for this module hangs from the branch containing the gufunc harness fix

Contributor

ovillellas commented Feb 20, 2013

At this point I think I've catered with all feedback in the comments: this means:

  1. Made sure that spacing in python functions followed PEP8 recommendations.
  2. Modified a bit documentation removing some "fails" on my side.
  3. When a NaN is generated I raise "invalid". This allows changing the behavior by using np.seterr. Updated doc accordingly.
  4. Fixed a build error when compiling for Python3 (detected by Travis)
  5. Filled the missing doctests.

I also found a bug in eig when using complex parameters and fixed. Also updated docstrings making sure that all cholesky related functions specify that they work for symmetrical/hermitian positive-definite matrices (in some places only symmetrical/hermitian was stated, which is not enough).

@ewmoore ewmoore and 2 others commented on an outdated diff Feb 21, 2013

numpy/core/src/umath/gufuncs_linalg.py
+ ----------
+ a : (<NDIMS>, N) array
+ Input array
+ b : (<NDIMS>, N) array
+ Input array
+
+ Returns
+ -------
+ inner : (<NDIM>) array
+ Inner product over the inner dimension.
+
+ Notes
+ -----
+ Numpy broadcasting rules apply when matching dimensions.
+
+ Implemented for types single and double. Numpy conversion rules apply.
@ewmoore

ewmoore Feb 21, 2013

Contributor

Why not also implement this for complex numbers? The inner product is well defined. Also a quick look seems to show that inner1d and innerwt are the only functions that don't work for complex arguments. It would be nice if they all could, since otherwise there is some asymmetry that seems unnecessary and potentially confusing.

@ovillellas

ovillellas Feb 21, 2013

Contributor

We didn't need them at the time, and there were questions about the proper way to go with this.
In fact, these functions came out as we were using them from the umath_tests, where inner1d and innerwt are implemented for float32 and float64. We wanted an accelerated version of inner1d implemented using BLAS sdot and ddot. We also made a version of innerwt to have both together.

We thought about adding a complex version. In BLAS there are two different versions of dot for complex numbers. [cz]dotu and [cz]dotc. Probably is dotu what would make sense, but as there was not a clear choice we just skipped this. If there is a consensus on which one to use that would be great (and quite easy to add). There is the option of supporting both with, which would involve a bit of extra effort. In any case, it was even less clear which would be the right option in the case of innerwt.

Really open to suggestions in this point :)

@charris

charris Mar 2, 2013

Owner

We really need both. A cdot might be a good addition to numpy.

@ovillellas

ovillellas Mar 15, 2013

Contributor

I am adding inner1d complex versions that will map to blas _dotu. I am adding also a tentatively named dotc1d that maps to _dotc. dotc1d and inner1d will map to the same function for real types. I will also expand innerwt to complex types using the regular complex product (not the conjugated one).

Right now I am testing it.

Contributor

ovillellas commented Mar 1, 2013

Fixed a problem with eig functions that were causing problems (tests failing causing Travis to fail)

Owner

charris commented Mar 2, 2013

The 1.8 release notes need a description off this new functionality.They are in doc/release.

I'm not sure whether this is the correct place to report, but I'm getting the following error:

import numpy as np
import numpy.core.gufuncs_linalg as gula
A = np.arange(2*2).reshape((2,2))
B = np.arange(2*1).reshape((2,1))
gula.matrix_multiply(A, B)
---------------------------------------------------------------------------
ValueError: On entry to DGEMM parameter number 10 had an illegal value
Owner

charris commented Apr 3, 2013

Commenting to ping the mail among other things, new commits don't do it.

Looking at this there are two big parts: changing lapack_lite and one adding gufunc. Would it be possible to break this into two PR's, one for each? I don' t have any objection to changing lapack_lite, in fact I'd like to see more of the existing functions exposed in numpy, but I'd like it to go in as a separate PRs.

I think the easy way to separate things (and squash the commits) is to check out a new branch, then checkout the modified files you want from this branch and commit.

Owner

charris commented Apr 3, 2013

Also, the new blas_lite and lapack-lite are in core/src/umath with only C exposure and that probably needs some discussion. Maybe it should be a library? Maybe we should move the originals in there also? Just saying it might be nice if everything was in one spot, python exposure of not.

Owner

charris commented Apr 4, 2013

I'm going to raise the library issue on the list.

@charris charris and 1 other commented on an outdated diff Apr 4, 2013

numpy/core/code_generators/cversions.txt
@@ -13,3 +13,5 @@
0x00000007 = e396ba3912dcf052eaee1b0b203a7724
# Version 8 Added interface to MapIterObject
0x00000008 = 17321775fc884de0b1eda478cd61c74b
+# Version 9
@charris

charris Apr 4, 2013

Owner

I don't think we are at 9 yet.

@ovillellas

ovillellas Apr 5, 2013

Contributor

In fact, I remember changing this to get around a build problem. I will revert that back, as it doesn't seem necessary

@charris charris commented on an outdated diff Apr 4, 2013

numpy/core/src/umath/gufuncs_linalg.py
+ 'multiply4_add', 'eig', 'eigvals', 'eigh', 'eigvalsh', 'solve',
+ 'svd', 'chosolve', 'poinv']
+
+
+import numpy.core._umath_linalg as _impl
+import numpy as np
+
+
+def inner1d(a, b, **kwargs):
+ """
+ Compute the dot product of vectors over the inner dimension, with
+ broadcasting.
+
+ Parameters
+ ----------
+ a : (<NDIMS>, N) array
@charris

charris Apr 4, 2013

Owner

The way to indicate shape is

a : array_like, shape (..., N)

I used array_like here instead of array as I assume nested lists and such will be converted. Is that so?
Needs fixing in the rest of the documentation.

@charris charris commented on the diff Apr 4, 2013

numpy/core/src/umath/gufuncs_linalg_contents.rst
@@ -0,0 +1,104 @@
+=======================
@charris

charris Apr 4, 2013

Owner

Where does this documentation go? I don't see it imported as part of a module anywhere or listed anywhere in doc/source/reference. I think you should put it somewhere in doc/source/reference, as undocumented features are ghostly creatures whose existence is in doubt.

@charris charris commented on the diff Apr 4, 2013

numpy/core/src/umath/gufuncs_linalg.py
+This module itself provides wrappers to kernels written as numpy
+generalized-ufuncs that perform the heavy-lifting of computation. The wrappers
+exist in order to provide a sane interface, like providing keyword arguments in
+line with the ones used by linalg, or just to automatically select the
+appropriate kernel depending on the parameters. All wrappers forward the
+keyword parameters to the underlying generalized ufunc (the kernel).
+
+The functions are intended to be used on arrays of matrices. Functions that in
+numpy.LinAlg would generate a LinAlgError (for example, inv on a non-invertible
+matrix) will just generate NaNs as result. When this happens, invalid floating
+point status will be set. Error handling can be configured for this cases using
+np.seterr.
+
+Additional functions some fused arithmetic, useful for efficient operation over
+"""
+
@charris

charris Apr 4, 2013

Owner

Could you update the new python files by adding

from __future__import division, absolute_import

at the top, but after the module docstring? It will help keep them compatible with the ongoing 2t03 project.

Owner

charris commented Apr 5, 2013

I'm about ready to put this in, it doesn't touch the existing files very deeply and the library issues can be hashed out later. The documentation needs some work, in particular you should put something in doc/source/reference, maybe the existing rst file with in entry in an appropriate index. An explanation of the new features should also be added to doc/release/1.8.0-notes.rst. To check that the documentation looks right you should also build it. Do cd into doc/sphinxext and install it, then in doc do make html, you can then view the result in firefox or the browser of your choice.

Owner

njsmith commented Apr 5, 2013

Now that the error reporting and such have been fixed (I think?), I'd like
to ask again what the reasons are that these functions can't generally
replace the ones in np.linalg? Are there any?

My worry is that we have a bad habit of just piling new stuff on top of old
with partial duplication of function and eventually ending up with tangled
deprecation cycles. And I don't really believe anyone will go back to clean
this up in a second pass if we merge a minimally usable version, so now is
the time to at least have the conversation.
On 5 Apr 2013 01:04, "Charles Harris" notifications@github.com wrote:

I'm about ready to put this in, it doesn't touch the existing files very
deeply and the library issues can be hashed out later. The documentation
needs some work, in particular you should put something in
doc/source/reference, maybe the existing rst file with in entry in an
appropriate index. An explanation of the new features should also be added
to doc/release/1.8.0-notes.rst. To check that the documentation looks
right you should also build it. Do cd into doc/sphinxext and install it,
then in doc do make html, you can then view the result in firefox or the
browser of your choice.


Reply to this email directly or view it on GitHubhttps://github.com/numpy/numpy/pull/2954#issuecomment-15931835
.

Contributor

ovillellas commented Apr 5, 2013

I will make the changes in the python module, and put the reference in doc/source/reference, as well as adding some release notes.

About blas-lite and lapack-lite, I agree that it would be nice to have a single version, probably with all the bindings for the new functions as well. However I think thats better done in a consolidation step. Right now the version of blas-lite/lapack-lite in umath is just statically liked to resolve the symbols in the c-module. It is only used if numpy is not configured with another blas/lapack provider (MKL/ATLAS/Accelerate on Mac...). Probably the right way would be to just have one lapack-lite library built with the symbols needed by both components.

About replacing the old linalg.py. In that sense the gufunc versions are not complete as to provide a full replacement for linalg.py. It could eventually if enough effort is thrown it, but right now I don't think it is possible. Apart from that, there are some functions that cannot be made to operate the same way that linalg versions do. For example the eig functions for real operands return real results in linalg if there all the eigenvalues are real. There is no way to do that in a gufunc in a satisfactory way, as a type for the result is needed before solving the problem, and because several problems could be paired together.

So the reasons for not being a replacement (right now) to linalg are completeness and backwards compatibility. The first problem will require some extra effort, the second one will need discussion.

Owner

njsmith commented Apr 5, 2013

It sounds like those issues (completeness and compatibility) are reasons why we can't just delete linalg.py and replace it with this module. But that doesn't mean they need to remain totally separate...

Can't we replace the existing linalg bindings with these for the functions that are in both (sticking to the old code for the functions that are in both)? And for compatibility issues, isn't this why the new code already has python-level wrappers for many of the new functions, so we could stick to that?

I just don't see any good reason to have two separate C-level wrappers for things like svd, nor any reason that it would cause problems for the existing np.linalg entrypoints to start transparently handling higher-dimensional arrays.

Contributor

ovillellas commented Apr 5, 2013

The only code that is really repeated is the blas-lite and lapack-lite source code that is just a f2c version of the stantard blas/lapack routines. Note that this is only used as a fallback when there is no other blas/lapack provider. As charis mentioned this code could be merged (it is mostly a grunt job involving regenerating the source code using f2p using the union of the functions in both modules, and modifying build scripts).

Other than that, there are some differences in the way both modules are implemented that make hard using the kernels in gufuncs in linalg transparently:

  1. linalg uses wrappers to the lapack functions directly, gufunc_linalg uses kernels written in C that call the lapack functions as part of the execution. That C code implements adapter code that in linalg is written in Python.
  2. Error handling behavior would be hard to reproduce. The checking of return values in Lapack functions in performed in C in gufunc-linalg. The condition is unavailable to the Python code. The fix to the issue with exceptions makes a "NaN encountered exception" if there is a NaN generated in any of the elements. On a Lapack error you will get NaNs and a "NaN exception" while in Linalg you will get a LinAlgException with a description of the Lapack error code. It is quite different and may not be possible to adapt.
  3. Linalg always uses double precision, independently of precision in input parameters. The gufunc_wrappers use single precision lapack functions for single precision operands (this behavior could be wrapped though).
  4. Linalg is probably nicer with matrices, although this could be wrapped as well.

It may be possible to implement some linalg functions that can not raise exceptions via gufunc-linalg. Maybe the differences in error handling are not important compared to the benefits of having the functions broadcast. I guess it is something to be discussed. But I also think that it would be a good idea to have the new kernels available for a while before modifying a module that has been working for a while, and consolidate it later after some additional usage.

Owner

njsmith commented Apr 5, 2013

linalg uses wrappers to the lapack functions directly, gufunc_linalg uses kernels written in C that call the lapack functions as part of the execution - obviously if linalg used the gufuncs as its backend, then this would change. I don't see how it matters how the old code that we'd be throwing away works.

Error handling behavior would be hard to reproduce - If the gufunc code has bad error handling, then this is just a bug in the gufunc code that should be fixed regardless. I suppose the best fix would be to enhance the general (g)ufunc machinery so that error reporting is possible -- this is a pretty fundamentally useful thing to have! -- but lacking that, as a for-now hack you could easily implement is to squirrel the Lapack error code away in a private TLS variable and check for it afterwards. (See PyThread_*_key in pythread.h for a TLS implementation that's always available. The documentation is in the Python source distribution in comments in Python/thread.c.)

Linalg always uses double precision, independently of precision in input parameters. The gufunc_wrappers use single precision lapack functions - This is exactly the kind of subtle incompatibility BS that makes me opposed to having duplicate APIs! Either linalg is wrong, or the gufunc_wrappers are wrong, and one of them should be fixed. Going forward, this kind of stuff will just crop up repeatedly, as people make enhancements to one or the other, or discover more subtle incompatibilities. Just say no to subtly incompatible APIs.

Linalg is probably nicer with matrices - not sure what "probably" means here :-) Surely linalg and the gufunc machinery both call __array_wrap__.

But I also think that it would be a good idea to have the new kernels available for a while before modifying a module that has been working for a while, and consolidate it later after some additional usage - No, I disagree here too. If we break linalg, we'll hear about it very quickly (during the RC cycle for anything even moderately bad), and it'll get fixed. If we stick these functions off in some ghetto, then they'll get minimal testing, and if gufuncs_linalg (or linalg) don't work for someone they'll just switch to using linalg (or gufuncs_linalg), and there will end up being people complaining how they were using one of them specifically because of some weird incompatible quirk (like the single-precision thing) and how dare we fix this, and it will be a nightmare to unravel. This kind of thing just makes more work for maintainers, and we don't have enough maintainers. Have the courage of your convictions! Working code is working code!

Contributor

ovillellas commented Apr 8, 2013

Error handling is different for a reason. You can solve many problems with a single call. That means that many can fail with different condition codes. Design was solve all the problems regardless of the error conditions of any single element and provide a way to detect the problem. If you have an array of say 10000 eigvalues problems in an array and the last one fails, you don't want to throw away 9999 valid results. And then, should it report the first error? or all the errors? what happens with the valid results?

IMO, In order to have "good" error handling there should be some machinery to handle per element error handling and correction.

About matrices... you are right, I guess the machinery wraps matrices properly (or at least if should). there is code in linalg to handle wrapping, but that's because it creates its own temporary buffers.

I just want to add: The idea of having this module was two-fold: convenience, as it allows broadcasting to work and performance, as it eliminates any need of looping in python code and removes some python glue code. The second point is performance.
In that context:

  • Marking errors instead of raising exceptions was selected as it allows not throwing away already computed results.
  • Using single precision to solve the single precision problems. They can be quite faster in the presence of optimized Lapack/Blas

The library was designed with that in mind, not to replace linalg. In that sense it adds value by adding a new tool, not by replacing an existing one. In fact, both complement each other, as an error in an entry executed via gufunc_linalg (that can be detected) could use code in linalg to diagnose the error, having the best of both worlds (speed when everything runs fine, better error diagnosis when something goes wrong.

Owner

charris commented Apr 8, 2013

In thinking about the library, Numpy tradition is to support a limited number of operations in double precision, while Scipy provides full access to LAPACK. In light of that, I think it would make sense to make the modifications to the ufunc machinery, select among the new c functions to provide basic functionality in Numpy, and move the rest to Scipy where the presence of the LAPACK library can be assumed. @rgommers @pv What do you think about that?

As to error handling, IMHO it has always been a bit inadequate in the ufunc loops and if this work leads to improvement so much the better. It might be worth bringing that up on the list for discussion.

I'd like to get much of this into Numpy without too much hassle. It is difficult to design and implement a perfect system at the first go, usually it works better to get something out there for people to work with and gather some feedback. It is a variant of the release early, release often mantra.

Owner

pv commented Apr 8, 2013

I think in the long run it's important to get rid of the duplication of work vs. scipy.linalg and numpy.linalg --- there's little point to have the two compete against each other in features (like this one), that's just wasted effort. numpy.linalg is there now, and we can't get rid of it, so I don't like the approach where Numpy contains "lite" versions of the routines --- there's no obvious upside, and it divides the number of people working on each by two (round down).

One possible path to success is to get the version of the functions in Numpy to feature parity with those in Scipy, and after that just drop the Scipy versions (or make them thin wrappers only so that the rubber hits the road inside Numpy). So if a way to do single-precision arithmetic is added to the routines Numpy, that's good. Adding new routines probably should go through the question: is it basic linear algebra?

Having the gufunc versions be completely separate from the numpy.linalg versions is IMHO not very good, as the functional difference is that one of them can operate on a "stacked" set of problems and one cannot. The error handling issue can likely be worked around by defining the generalization from N=1 upwards reasonably (e.g. raise if everything fails, or just drop the concept of LinalgErrors as any code relying on them is likely numerically unstable) --- so I don't see fundamental problems with, say, numpy.linalg.solve becoming able to solve a stacked set of equation systems.

As far as this PR is concerned --- the feature itself is certainly something I've been wishing for. However, I don't think the API is yet in its final form. My take would be here: +1 for getting it in, but integrating it into numpy.linalg should be considered seriously. So, maybe merge and mark the feature as experimental, and we'll work on from that point?

As to the question whether something should go to Scipy: I think this decision should be made on the basis on which functions are "basic" enough to put in numpy.linalg. There's in principle nothing preventing gufuncs in Scipy, and the f2py based wrapping of BLAS/LAPACK currently used in Scipy is not an optimal way to go. I'd be happy to see Scipy linalg functions gufuncized.

(I think we can mostly ignore the LAPACK/BLAS issue --- in any serious environment Numpy is linked against a real LAPACK+BLAS combination, and we can always put more stuff to lapack_lite.c if necessary.)

Owner

charris commented Apr 8, 2013

Thanks Pauli. So if this is merged now the remaining tasks would be

  1. Combine the *_lite libraries.
  2. Expose more of the library contents in linalg.
  3. Use gufuncs in linalg.py where appropriate.
  4. Unify scipy/numpy linalg with the goal of moving all to numpy
  5. Work on the gufunc API and error handling.

My own priority would be 1), which is pretty much a requirement before proceeding with 2..4.

@ovillellas Oscar, I still don't see where the documentation goes. What is your long term interest improving numpy linalg?

Owner

rgommers commented Apr 8, 2013

@pv doesn't the history of this separation between numpy.linalg and scipy.linalg have to do with Fortran code? So are you proposing to drop only the scipy versions where there is an existing duplicate C version in numpy now? Therefore scipy.linalg will always contain things that won't be available in numpy.linalg, correct?

You're also talking about "basic enough". So we have a basic/non-basic and a C/Fortran divide. You're right that the duplication of work now is a waste. It would help me though to see an overview at some point of what goes where for the scenarios of (a) only dropping duplicates, and (b) gufuncizing everything.

Contributor

ovillellas commented Apr 8, 2013

My Idea was adding an entry to the sphinx documentation for the module, with bits from the gufuncs_linalg_contents.rst file. The idea would be remove that .rst in favor of the sphinx source. Maybe the more technical details could be bundled as a README that is more suited for development.

From discussions here I guess that in the documentation it should go in a way that makes clear the differences with linalg, and not to sell it as a replacement (that it is not right now).

As soon as I submit that I will ping this thread. I do think that in any case getting the doc write is quite important to avoid problems. I was placing a link to the new entry in the reference, inside the routines section, under the entry for linalg module. I was thinking naming it something in the lines of "Generalized ufuncs for Linear Algebra", but I am open on how to name that and where to place it.

WRT the question about my long term interest, we developed this as part of a project and wanted to share it as others may find use for it. I would help to better integrate the code basis, but will be limited on how much time I can put into it, depending on my workload on other projects.

I just wanted to add that while developing the module there were plenty of ideas around the gufunc interface. One was related to the error handling mentioned here, were exceptions didn't look like a good fit for vector executions of problems (as an exception kind of aborts the whole execution even if a single element fails). But there were other points raised (and that eventually resulted in the creation of the wrappers instead of it being just a collection of gufuncs in a C module):

  • The way signatures work constrains quite a bit what can be done in gufuncs, For example, in svd the signature should look like (m,n)->(min(m,n)). Of course, that is not supported. Probably is not worth supporting complex operations inside the signature, but having a signature like '(m,n)->(p)' and a function that is able to resolve p based on m and n could help. This could be part of an improved interface (having the ability to resolve and/or validate signatures). This was worked around by having two gufuncs: (m,n)->(m) and (m,n)->(n) and have the wrapper choose which one to use. But there may be cases where such an work around may not be feasible.
  • It would be useful a way to pass uniform parameters to gufuncs, That is, parameters that are shared among all the elements.
  • It would also be nice to have setup and tear-down functions for the whole execution. AFAIK, the machinery tries t o minimize the number of calls to the gufunc kernel, but it is not always possible. In the current machinery initialization and tear-down of working buffers needs to be done per-kernel call, which can happen several times for a single call. The buffers could be reused if such setup/tear-down functions existed. However, this is not that critical as it won't happen in the inner-most loop anyways (and the machinery tries to minimize the number of calls as well).
Owner

pv commented Apr 8, 2013

@rgommers: the historical reason is probably to keep lapack/blas_lite.c reasonably small. The interface to Lapack goes anyway through the Fortran ABI, and numpy.linalg functions indeed link against the system LAPACK/BLAS when available. Wrapping all of LAPACK (well, at least the SUBROUTINEs) should be possible in C.

Also yes, I was only proposing to drop the already existing duplicate functions --- this seems like the path of least resistance. There's additional stuff in Scipy like solve_lyapunov etc. which I don't believe belong in Numpy. Perhaps there's some stuff that should be in Numpy but isn't, can't say immediately.

EDIT: forgot, expokit is eventually coming to scipy.linalg, so there will indeed be more Fortran code there in the future.

Owner

njsmith commented Apr 8, 2013

I agree with Pauli that the numpy.linalg/scipy.linalg distinction is not
anything we want to put resources into preserving, given that it's only
reason for existing is to make developers lives somewhat easier. I can't
imagine any users actually benefit from having two sets of near-identical
APIs...

Regarding single vs. double precision, I also think that if single
precision is useful, then it should be available in numpy.linalg -- so
maybe we ought to bite the bullet and transition the numpy.linalg routines
to respecting input precision by default? (A nice thing about reusing the
ufunc machinery would be that if anyone wants to use double-precision
intermediates for single-precision inputs, there are already mechanisms for
this via dtype=, loop=, etc. Anything else requires weird and inconsistent
kluges.)

Regarding error handling: it would indeed be excellent to have richer error
handling for gufuncs than a simple exception/no-exception. There seem like
a lot of possibilities for how this should actually work, though, that seem
like they need to be thought about a little -- do LAPACK routines ever
leave useful partial results behind when they error out? is the actual
error code useful? Should we re-use the seterr() mechanism, or just raise
an exception with some extra information appended (so people who want to
catch the exception can extract the partial results from it), or something?
If we use seterr should we really re-use the 'invalid' error code for this
(as opposed to defining a new one), given that it's supposed to indicate
some very specific conditions defined in the IEEE754 spec? My guess is that
these are actually probably not too hard to figure out if we think about it
but we'd have to do the thinking.

I take Chuck's point about not trying to be perfect on the first go, but if
we just merge this as an "experimental" feature in a new parallel namespace
somewhere, then here's what I predict will happen:

  • Initial itch scratched, we'll all go on to other stuff and this'll pretty
    much sit modulo the occasional patch for a few years
  • During that time, we'll get a steady stream of users confused about what
    the difference between the two (or three, including scipy) linalg routines,
    and frustrated trying to figure out which one they actually need
  • There'll be random bugs and inconsistencies; patches will either be
    applied to just one of the trees, or we'll have to go to extra effort to
    apply bug-fixes repeatedly. (E.g., think of the ongoing qr changes.)
  • Eventually, this will get annoying enough that someone will take the
    effort to clean things up, but this will be extra difficult because now
    there will be code in the wild that depends on the unique quirks of each of
    these interfaces, so the final unified code will have to go to extra work
    to emulate the bugs in each of these old interfaces during some deprecation
    period.

Release early, release often is awesome, but we still have to think about
quality sooner or later, and sticking an "experimental" label on things
doesn't magically mean we don't have to worry about backwards
compatibility. IMO when there are API-level things that we know will need
fixing later, and where a bit of thought actually will help, then we
should fix them before merging instead of cutting off the conversation. (I
say "API-level", because those are the things that end up creating
compatibility entanglements. The BLAS/LAPACK duplication isn't nice, but at
least ignoring it for now won't create any new problems down the road, so I
don't care as much.)

This will all be much nicer if we can (1) make sure that this code is
actually used by the standard interfaces that everyone knows and uses,
which will mean that it gets way more loves and testing, (2) avoid
publically exposing APIs that we haven't even thought through. (Putting
stuff into the wild is a good experiment, but like every experiment you
really have to think out all the parts you can before-hand if you want to
learn anything useful.)

To show how close we are, here (I believe) is a complete, strictly
compatible, more featureful implemention of np.linalg.solve, assuming the
gufunc code as it currently exists:

def err(_args):
raise LinAlgError
def solve(a, b, *args, *_kwargs):
kwargs.setdefault("dtype", np.float64)
with np.errstate(invalid="call", call=err):
return solve_ufunc(a, b, _args, *_kwargs)

This is truly trivial, and we could easily phase out the compatibility
hacks over time (via deprecation periods etc.).

On Mon, Apr 8, 2013 at 4:28 PM, Pauli Virtanen notifications@github.comwrote:

I think in the long run it's important to get rid of the duplication of
work vs. scipy.linalg and numpy.linalg --- there's little point to have the
two compete against each other in features (like this one), that's just
wasted effort. numpy.linalg is there now, and we can't get rid of it, so
I don't like the approach where Numpy contains "lite" versions of the
routines --- there's no obvious upside, and it divides the number of people
working on each by two (round down).

One possible path to success is to get the version of the functions in
Numpy to feature parity with those in Scipy, and after that just drop the
Scipy versions (or make them thin wrappers only so that the rubber hits the
road inside Numpy). So if a way to do single-precision arithmetic is added
to the routines Numpy, that's good. Adding new routines probably should go
through the question: is it basic linear algebra?

Having the gufunc versions be completely separate from the numpy.linalg
versions is IMHO not very good, as the functional difference is that one of
them can operate on a "stacked" set of problems and one cannot. The error
handling issue can likely be worked around by defining the generalization
from N=1 upwards reasonably (e.g. raise if everything fails, or just drop
the concept of LinalgErrors as any code relying on them is likely
numerically unstable) --- so I don't see fundamental problems with, say,
numpy.linalg.solve becoming able to solve a stacked set of equation
systems.

As far as this PR is concerned --- the feature itself is certainly
something I've been wishing for. However, I don't think the API is yet in
its final form. My take would be here: +1 for getting it in, but
integrating it into numpy.linalg should be considered seriously. So, maybe
merge and mark the feature as experimental, and we'll work on from that
point?

As to the question whether something should go to Scipy: I think this
decision should be made on the basis on which functions are "basic" enough
to put in numpy.linalg. There's in principle nothing preventing gufuncs in
Scipy, and the f2py based wrapping of BLAS/LAPACK currently used in Scipy
is not an optimal way to go.

(I think we can mostly ignore the LAPACK/BLAS issue --- in any serious
environment Numpy is linked against a real LAPACK+BLAS combination, and we
can always put more stuff to lapack_lite.c if necessary.)


Reply to this email directly or view it on GitHubhttps://github.com/numpy/numpy/pull/2954#issuecomment-16057166
.

Owner

pv commented Apr 8, 2013

@njsmith: agreed on the API stability --- your prophecy rings true. Getting as much as possible into numpy.linalg would go with the idea that there should be one --- and preferably only one --- obvious way to do things.

To keep the ball rolling: gufungs_linalg.py -> _gufuncs_linalg.py, and just merge without user-visible documentation. Then try to get some stuff into numpy.linalg, or argue that it should be advertised as an actual experimental API?

Owner

pv commented Apr 8, 2013

So, something like this seems possible: pv/numpy@a344d47...pv:linalg-gu
I think this fruit is hanging quite low. (Of course I picked the easiest functions first :)

Here, I punted on the float32 vs. float64 thing, and we may want to switch raising warnings instead of LinAlgErrors (at least for inversion), but things seem robust. Test suite passes whatever that means. If we want to improve the error handling mechanism, e.g. raising custom error types, then I think an extension of the current errobj mechanism will be a soft point for pushing more functionality in.

Overall, it seems to me that this functionality will be very easy to discover and understand if it's put in numpy.linalg.

Contributor

ovillellas commented Apr 9, 2013

@pv: Keep in mind that gufuncs_linalg.py is in itself a wrapper of _umath_linalg, the C extension module that implements the functionality. In many cases it is a straight call to the underlying C function, but I decided to add it as a full function so that docstrings for all functionality was in a single python file.

So if the functionality really ends into linalg, maybe gufuncs_linalg should just go and place the wrapping it provides in linalg itself. Or just keep it for the different exception semantics.

Another option would be allowing changing the error behavior by either a keyword parameter in the functions or by setting some variable in the module...
linalg.on_error='raise' (default, but could be set to 'ignore')

Owner

pv commented Apr 9, 2013

@ovillellas: yes, I think gufunc_linalg.py shouldn't be kept around when the same features are in numpy.linalg. I think it is not too difficult to collapse everything into numpy.linalg even when maintaining backward compatibility.

How the error handling should exactly be done can be discussed. One option indeed is to add global state, e.g. np.linalg.seterr which works similarly to np.seterr, or control it via a keyword argument. Or, extend np.seterr with a linalg keyword. Or, specify error handling via a keyword argument. Or, a combination of the above.

This is perhaps partly orthogonal to the generalized ufunc implementation, as also for single matrices it can be useful to signal errors with NaNs rather than LinAlgerrors. One possibility is to first merge the functionality into numpy.linalg raising LinAlgErrors, and then think how to allow users to e.g. convert errors to warnings, or ignore them altogether.

Owner

pv commented Apr 9, 2013

@charris, @njsmith: what do you think? Merge now, maybe cherry-picking 1a482c7 and df6fe13, and then do integration with numpy.linalg, sorting out how to get rid of LinAlgErrors and casting to doubles in separate PRs?

I have _umath_linalg based backwards-compatible (or so I hope :) versions of the common routines in the branch.

Owner

charris commented Apr 9, 2013

@pv I'm inclined to merge early, but I think it would be good to open a task, or list of tasks, for what remains. AFAICT, that is to unify the *_lite libraries, decide on and implement the error handling, and use the new routines in linalg as far as possible. It might be a good idea to sort the error handling before this goes in. However you want to get your work in works for me, either cherry-picking or a follow on PR would be fine.

Owner

njsmith commented Apr 10, 2013

@pv: Is it an option now to just merge your branch and this together? Would
that leave us with a working, compatible np.linalg? From a quick skim it
looks like you've just sorted out some small bug-fixes and have kept
np.linalg compatibility -- except I'm not sure about the single/double
issue, but that shouldn't be hard to sort out.

IIUC this would leave master in a working state where the core parts of
this code were actually in use, but adding a ton of new exposed API surface
area. There would still definitely be cleanups we wanted to do -- to code
layout, error handling, exposing more capabilities, etc. -- but I think at
that point they'd all be stuff we could do incrementally, without breaking
API in the process, and without risking this code getting lost.

On Tue, Apr 9, 2013 at 10:58 PM, Charles Harris notifications@github.comwrote:

@pv https://github.com/pv I'm inclined to merge early, but I think it
would be good to open a task, or list of tasks, for what remains. AFAICT,
that is to unify the *_lite libraries, decide on and implement the error
handling, and use the new routines in linalg as far as possible. It might
be a good idea to sort the error handling before this goes in. However you
want to get your work in works for me, either cherry-picking or a follow on
PR would be fine.


Reply to this email directly or view it on GitHubhttps://github.com/numpy/numpy/pull/2954#issuecomment-16143061
.

Owner

pv commented Apr 10, 2013

@njsmith: Yes, it was just replacing the core parts in numpy.linalg with the umath_linalg routines where applicable, keeping the old semantics for ndim <= 2, so nothing fancy. (The previous behavior with single vs. double issue is handled by converting single->double before processing and then back double->single.)

As far as I know, the backward compatibility is preserved as long as it's about arrays with ndim <= 2. So we could go at it incrementally. The test suite passes, but it should still be tested against some Numpy-using packages, though, didn't do that yet. Some 3rd party code might in principle depend on errors being raised for ndim > 2 input, but I think we can just go and break this.

Owner

njsmith commented Apr 10, 2013

The general rule is that we don't worry about error -> non-error
transitions when handling deprecations, yeah.

So maybe the best approach is for Pauli to go ahead and submit his branch
as a new PR (which will automatically include everything in this one,
because of how git history tracking works)? Does that make sense to
everyone else? @charris, @ ovillellas?

@njsmith https://github.com/njsmith: Yes, it was just replacing the core
parts in numpy.linalg with the umath_linalg routines where applicable,
keeping the old semantics. (The previous behavior with single vs. double
issue is handled by converting single->double before processing and then
back double->single.)

As far as I know, the backward compatibility is preserved as long as it's
about arrays with ndim <= 2. So we could go at it incrementally. It should
still be tested against some Numpy-using packages, though, didn't do that
yet. Some 3rd party code might in principle depend on errors being raised
for ndim > 2 input, but I think we can just go and break this.


Reply to this email directly or view it on
GitHubhttps://github.com/numpy/numpy/pull/2954#issuecomment-16183273
.

Contributor

ovillellas commented Apr 10, 2013

@njsmith Makes sense, I guess the discussion should be linked somehow, just adding a pointer to the new pull request will suffice.

Pauli's changes just modify linalg itself, right?

Just a note: in gufuncs_linalg there is some extra function not present in linalg, so it may be a good idea to add that as well. I am thinking of chosolve (solve using cholesky decomposition for Hermitian, positive-definite matrices).

Owner

njsmith commented Apr 10, 2013

@ovillellas: Pauli's changes are here, you can look for yourself :-)
pv/numpy@a344d47...pv:linalg-gu

Looks like they also add an _ to the beginning of the name of
gufunc_linalg.py, and add some extra FP flag clearing to umath_linalg.

chosolve sounds like an excellent addition indeed...

On Wed, Apr 10, 2013 at 5:41 PM, Óscar Villellas Guillén <
notifications@github.com> wrote:

@njsmith https://github.com/njsmith Makes sense, I guess the discussion
should be linked somehow, just adding a pointer to the new pull request
will suffice.

Pauli's changes just modify linalg itself, right?

Just a note: in gufuncs_linalg there is some extra function not present in
linalg, so it may be a good idea to add that as well. I am thinking of
chosolve (solve using cholesky decomposition for Hermitian,
positive-definite matrices).


Reply to this email directly or view it on GitHubhttps://github.com/numpy/numpy/pull/2954#issuecomment-16186285
.

Owner

pv commented Apr 10, 2013

The PR is here: gh-3220 and a task issue is here: gh-3217

Yes, there are several new functions left that would be useful to add into numpy.linalg, but I left that for later.
At least for chosolve some though needs to be paid to the corresponding Scipy function signatures...

Owner

njsmith commented Apr 11, 2013

Closing this to keep the PR list tidy (or tidier...).

njsmith closed this Apr 11, 2013

@larsmans larsmans commented on the diff Sep 30, 2013

numpy/core/src/umath/gufuncs_linalg.py
+on the operands.
+
+This module itself provides wrappers to kernels written as numpy
+generalized-ufuncs that perform the heavy-lifting of computation. The wrappers
+exist in order to provide a sane interface, like providing keyword arguments in
+line with the ones used by linalg, or just to automatically select the
+appropriate kernel depending on the parameters. All wrappers forward the
+keyword parameters to the underlying generalized ufunc (the kernel).
+
+The functions are intended to be used on arrays of matrices. Functions that in
+numpy.LinAlg would generate a LinAlgError (for example, inv on a non-invertible
+matrix) will just generate NaNs as result. When this happens, invalid floating
+point status will be set. Error handling can be configured for this cases using
+np.seterr.
+
+Additional functions some fused arithmetic, useful for efficient operation over
@larsmans

larsmans Sep 30, 2013

Contributor

Unfinished sentence (also missing verb, typo in "additional")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment