ENH: re-enabled old sv lapack routine as defaults

[zerothi]

This fixes the speed regression in #7847. Since the default
change of sv to svx the solve routine suffered a huge performance
penalty for anything but low order NRHS.

This commit fixes that issue by enabling the use of either 1) the
svx or 2) the sv routines. With a default on the latter.

However, the main funtionality by using svx was the easy check
of the condition number to assert a non-singular matrix.
This commit adds a call to the con routines to extract the
appropriate condition number. This force the addition of:

• lamch (machine precision extraction)
• gecon (condition number calculation from LU factorization)
• lange (1-norm of matrix)

This commit adds to arguments (and deprecates one) by:

• refine, bool, decides whether the svx (true) or sv (false)
  routines should be used.
• form = 'none', 'trans', 'conj' which refers to the form
  of the solution step. This is only valid if refine=True
  A ValueError will be issued if not none and refine=False.
  This keyword deprecates the "transposed" boolean value which
  only allowed 'T' (real) and 'C' (complex) but not 'T' for
  complex.

A fix for the deprecated sym_pos keyword which wasn't used in the previous
version.

A couple of additional tests have been added, mainly to check the arguments.

[zerothi]

Lets see if it passes tests through ci, it does so locally.

[pv]

Please make, as a first step, a PR that does not add any new options but only reverts to the gesv+gecon approach. I don't think we want to add them to 1.0.0 at the least --- it's already branched and only bugfixes (such as the one concerning the performance) should go in.

After that, we can discuss about the possibility to have the added options.

[ilayn]

Also there are symmetric but not hermitian complex matrices and "T" is used for these matrices for right inversion. Thus zsysv and zhesv solves different problems.

[zerothi]

@ilayn for the gesvx routines one can pass TRANS='T' or TRANS='C', for a^T x = b and a^H x = b, respectively. Both TRANS/='N' are only for the gesvx routines. Anything, but 'N' are not accepted for he/sy/po routines.

[zerothi]

@pv But the change to the other routines has a regression on the precision of the algorithm. The only way to circumvent this is by using the x routines of LAPACK.
Is this wanted?

[pv]

The *sv routines were used for many years without problems, so I doubt the precision issue is very important for most users.

[pv]

The issue is that there is not much time to hash out the design and any other issues in time for the release, so I'm -1 on adding new features to 1.0.0 at this stage.

[zerothi]

Ok, I will amend, this forces me to remove the transposed keyword (issue a deprecation warning) as it is only available for the gesvx routines.

[rgommers]

@zerothi there's nothing wrong with this PR for master, so how about leaving the desired behavior in here, and sending a separate PR against 1.0.x that simply reverts to the pre-0.19.0 situation?

[zerothi]

ok, will do. Will then just add a deprecation warning for the the transposed keyword in this pr to make it follow the PR i will create for 1.0.x.

[pv]

Note that it's not possible to remove the `transposed` keyword in 1.0.0 since we'd have to deprecate it first. It's possible to implement it via calling getrf+getrs instead of gesv. That's actually what gesv does, http://www.netlib.no/netlib/lapack/double/dgesv.f

[zerothi]

As I said, I won't remove it, simply deprecate it?

[zerothi]

What I can do is:

for master enable getrf + gesv with my newly changed form keyword

1. This current PR with some minor changes to allow form for assume_a = 'gen', regardless of used routine (gesv vs. gesvx).
2. Deprecate transposed keyword

for 1.0.x, do not allow gesv and transposed keyword. Its use only allowed T for real, and C for complex, it didn't allow T for complex cases. So if you want the same functionality you are imposing certain constraints on the user, i.e. the user may expect to get T for a complex solution, whereas it actually does C?
My point being, transposed=True isn't a ^T x = b when a is complex in the current solve.

Or would you rather that transposed in 1.0.x behaves as in the old version, which in my opinion will just be a serious regression when they go to 1.X>0.Y?

[pv]

"Luckily", there's a bug for complex so there are no backward compatibility issues for what transposed means in that case: ``` >> scipy.linalg.solve(np.eye(3).astype(complex), np.ones(3), transposed=True) ** On entry to ZGESVXUMEpsilonFull parameter number 2 had an illegal value ValueError: LAPACK reported an illegal value in 2-th argument. ``` For 1.0.0 probably best to just explicitly raise an error for the complex transposed case. For 1.1.x we might want to reuse keyword argument name, `transposed={True, 'C', 'T', 'N'}`, so that it's not necessary to deprecate the argumetn name. The behavior in the complex case for transposed=True can remain an error since it's ambiguous which one was meant --- or it should be made to be 'C' to remain consistent with assume_a='her'.

[zerothi]

@pv i think the transposed keyword was passing 'H' in the old scipy, whereas it only accepts 'C'. Hence the argument error. So yes, a bug in the previous solve for complex cases.

As for 1.1. Since the transposed keyword is ambiguous, wouldn't it be better to rename it to form? I mean the form of the equations to be solved is non-ambiguous whereas transposed is not so much, could be a^T, or x^T.

I will hardly think many people are aware of its existence (it just entered in 0.19) and we can retain it in 1.1 as a backfall (same as sym_pos keyword).

I much prefer fixing non-ideal keyword names, rather than continuing with ambiguous names?

So, for master, I know what to do, for 1.0.x, please suggest your preferred way of route, still transposed={False, 'N', 'T', 'C'}, or as I suggested: form={'none', 'trans', 'conj'}?

[ilayn]

I'm having hard time following the discussion but I think I'm a bit closer to what @pv mentioned about reverting it to old version (if it is really the problem, since we also didn't have any complaints regarding the speed until now). Because there is a time pressure for 1.0 or not touching it and performing a full surgery until 1.0.1 or 1.1. Yesterday, I've done some tests and indeed the time is mostly spent on refinement steps.

Keyword names are of course open for improvement or renaming. But transpose issue is simple: complex matrices can also have T keyword for complex generic matrices and if sym is given for complex matrices it shouldn't switch to hermitian solution. Otherwise it is not correct.

I didn't test line by line but just mentioning it if that's the case. But form is not related to the thing it represents or maybe I don't understand what form refers to. It's a switch between AX=B and A^(T)X=B so it's a property of the coefficient matrix. Also LAPACK naming is similar with trans keyword.

[ilayn]

I'll check about the bug as soon as I have some time.

[zerothi]

@ilayn yes, I agree to the TRANS name, however the lapack documentation specifies that the equation has the form (hence the name) which I think justifies the name, short and concise.
If there are other, better names, I am very much open to suggestions ;)

[zerothi]

I have now squashed some testing commits into two commits and here are some final remarks that needs either approval or change, but at least requires some comments from you:

1. When reverting to the gesv routines certain tests was failing. One was fixed in c7370f3 which I think can easily be justified from the documentation of the linprog routine, i.e. it automatically backfalls to the non-symmetric case. This wasn't seen when using the gesvx routines because the sym_pos keyword was affected by the bug as mentioned above. I can't see this should be a problem.

2. Another error arises in the TestExpM.test_triangularity_perturbation (scipy/sparse/linalg/tests/test_matfuncs.py)
   Here a RuntimeWarning is issued due to the condition number. In all older versions of scipy, (pre gesvx) this could not have been issued due to no checks. For the 0.19.X versions the gesvx routines actually passes these tests and thus it was a hidden failure. I am not sure what you think should be the path forward? I am not experienced enough to see if this is due to a corner case of the created test, or whether the _solve_P_Q should fall-back on the gesvx routines?

Could you comment on the above things?
Sadly, when backporting this to 1.0.X the second option will be a problem because you intent to remove gesvx. How should we go about that?

[pv]

For the expm case, the bad conditioning possibly doesn't matter for the algorithm, so the warning can be ignored. This is easiest to do in the test itself, by ``` from scipy._lib._numpy_compat import suppress_warnings ... def test_triangularity_perturbation(self): ... with suppress_warnings() as sup: sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll-conditioned.*") A_expm_logm_perturbed = expm(A_logm_perturbed) ... ```

[zerothi]

Ok, I can add that, but what about 1.0.X?

[pv]

Also, the expm test was written well before the change to gesvx, so I believe the test works as intended regardless of the bad conditioning.

[zerothi]

Ah, yes, of course. I see. I will do that!

[zerothi]

@pv @rgommers @ilayn ready for review

[ilayn]

@zerothi I think after the cherrypicking the commits this is also now redundant right?

[zerothi]

Agreed, this PR was trying to retain both vx and v routines. Since the discussion landed on not needed I will close, thanks!