pytensor icon indicating copy to clipboard operation
pytensor copied to clipboard

Directly use scipy lapack functions instead of scipy wrappers

Open jessegrabowski opened this issue 10 months ago • 6 comments

Description

We take small performance hits by using scipy.linalg functions rather than directly calling lapack functions in our perform methods. These wrappers include a lot of data validation checks that we don't need, and sometimes end up copying things that we don't want to copy.

Instead of using these wrappers, we can use scipy.linalg.get_lapack_funcs to get and use the relevant routine. For example, calling potrs instead of using linalg.cholesky.

Which routine to call is easy to figure out by looking at the scipy source code and checking which function the wrapper calls, and with which arguments.

Here's a list of Op that can be converted, and the LAPACK routines they call:

  • [x] Cholesky (potrf)
  • [ ] CholeskySolve (potrs)
  • [ ] LUFactor (getrf)
  • [ ] LUSolve (getrs)
  • [ ] SolveTriangular (trtrs)
  • Solve (depends on the assume_a argument. They all do lu_factor + lu_solve, but the exact routines called depends on the structure of the A matrix)
    • [ ] "gen" - getrf + getrs
    • [ ] "sym" - sysv
    • [ ] "pos" - posv
    • [ ] "tridiagonal" - gttrf + gttrs
    • [ ] "banded" - gbsv
  • [ ] Eigvalsh (complicated -- depends on input arguments. Need to spend some time looking at the scipy wrapper)
  • [ ] SolveContinuousLyapunov - trsyl

Routines in tensor.nlinalg might also benefit, but it would require some testing. For example, GESVD and GESDD are LAPACK routines for SVD, which might be better than np.linalg.svd (and related functions, like np.linalg.pinv). As I type this, I can't say if calling these routines directly would be better or worse than calling the numpy function. The numpy code is somewhat harder to peer into than scipy.

jessegrabowski avatar Jun 12 '25 23:06 jessegrabowski

I would like to work on this.

Would it be necessary/relevant to replace the wrappers outside of perform methods? (e.g. lu_factor.py, test suites )

Fyrebright avatar Jun 18 '25 00:06 Fyrebright

No, this issue only concerns the perform methods.

I suggest you open a PR to just handle one Op, and we can iterate from there. I think Cholesky will be a relatively simple one to start with.

jessegrabowski avatar Jun 18 '25 19:06 jessegrabowski

I am unsure if I've missed some assumptions that would eliminate more of the checks that are in the scipy wrapper. There seems to be a small improvement.

The check_finite behavior/default was inconsistent with the docstring so I changed it to false, but maybe the documentation should have been changed instead?

Fyrebright avatar Jun 19 '25 22:06 Fyrebright

I updated this issue to have a check-list of functions that could benefit from having their wrappers stripped away.

Benchmarking on the routines that would replace tensor.nlinalg functions would be appreciated.

jessegrabowski avatar Jun 23 '25 09:06 jessegrabowski

Since the changes were very similar to Cholesky, I thought it would be more convenient for reviewing to put CholeskySolve, LUFactor, SolveTriangular in one PR. If not I can make separate ones for each.

Fyrebright avatar Sep 02 '25 16:09 Fyrebright

The existing lu_solve/_lu_solve functions don't call Scipy and just use slinalg.solve_triangular. Should this still be changed to use getrs if solve_triangular is no longer using Scipy?

Fyrebright avatar Sep 02 '25 17:09 Fyrebright