Re: Another problem with matrix inversion

[Bill Allombert <Bill.Allombert@math.u-bordeaux1.fr>]

On Mon, May 25, 2009 at 10:22:39AM +0200, Lorenz Minder wrote:
> Hi,
> 
> BA:
> > On Thu, May 21, 2009 at 04:43:46AM +0200, Lorenz Minder wrote:
> > > > Sorry for the rushed our reply, but maybe you would prefer it to
> > > > a long delayed one.
> > >
> > > I absolutely do, yes.  Thanks for taking the time.
> > 
> > But I am afraid this lead to misunderstanding. Sorry about that.
> 
> Maybe, and with my rushed reply on your reply, it probably has irritated
> both of us to some degree.  In any case, I apologize form my harsh tone
> in that last mail.
> 
> I still think your "rushed reply" was a good thing, after all it helped
> clear out the non-controversial issues.
> 
> > > I'll look in the test suite, not sure how this works.  (Documentation
> > > is noted, let's talk about this later.)
> > 
> > Test-suite is merely a matter of adding a file src/test/in/foo or
> > expanding an existing file that is not used by 'make bench'.
> 
> Ok.  Is in/linear used by `make bench'?  How do I know?

%make bench|grep linear
* Testing linear        for gp-sta..TIME=20     for gp-dyn..TIME=20

so yes.

> As a tangential issue, when I ran the tests a couple of months ago I
> got a number of failures (Linux on PPC, GMP kernel IIRC), all of which
> were stack overflows.  Should I redo this (with unpatched source) and
> send you guys the results?  (At the time I didn't bother because it
> didn't look like something as serious as wrong output.)

Please do, however I have set up a buildlog on a PlayStation 3 which
does not show any problem currently.

> > > I strongly suggest to go back to the way I did it originally, namely
> > > simply replace FpM_gauss() and Flm_gauss().  Let me quote from the
> > 
> > I have an alternative:
> > 1) We add function FpM_invsafe/Flm_invsafe that would do that same for
> > matrix inversion that Fp_invsafe do for integers. Maybe even
> > FpM_gausssafe if needed.  This is a minimal change, that use a
> > well-defined interface.
> 
> I'm not too sure I get this right.  Fp_invsafe() gives an inverse if it
> exists, so FpM_invsafe would do the same (i.e., it would behave like
> ZnM_gauss(*, NULL, modulus))?  Or do you mean that FpM_invsafe does what
> FpM_gauss() does except err out when p is detected to be nonprime?
> Assuming that it then returns the failing residue class that would be
> useful, of course.

Well the general idea is to rename the current FpM_gauss to FpM_gausssafe,
change it to return the residue (t_INT) if it fails, or the matrix otherwise,
and have a wrapper FpM_gauss that do

GEN
FpM_gauss(GEN a, GEN b, GEN p)
{
  GEN z = FpM_gauss_safe(a,b,p);
  if (typ(z)==t_INT)
    pari_err(invmoder, gmodulo(z,p));
  return z;
}

> I think though that FpM_gausssafe() and Flm_gausssafe() are more
> important than the *_invsafe()-variants.  (Also *_invsafe() is just
> trivial if we have *_gaussafe().)

Agreed.

> > I can do it if you want.
> > 
> > 2) We add a function that take a factored modulus and use chinese
> > remainder and Hensel lifting to use FpM_invsafe to inverse a matrix.
> > Hensel lifting alleviate the need to change FpM_gauss to look for an
> > invertible pivot. This is similar to how gen_Shanks_sqrtn works.
> 
> Ok, so I guess it's the 2nd possibility above, right?
> 
> > 3) We add a function that take an unfactored modulus, call the
> > function above as if the modulus is prime, and if not, refine the
> > factorisation and iterate.
> 
> This sounds reasonable.  The only cost compared to the current solution
> seems to be that Gaussian elimination is restarted from scratch if the
> modulus is proved non-prime.  If that is the case, I think that would
> work just fine. (I don't worry too much about this added cost.)

Well, maybe step 2) could use step 3) for each assumed prime factor to
avoid this issue.

> I'm still not sure how your *safe() proposal would fix the problem with
> the proliferation of new Zn*-routines that the previous approach would
> cause.  Can you elaborate on that?

The idea is that ZnM_inv will defer all linear algebra to FpM_gausssafe,
so there will not need a lot of ZnM function anymore. We will still
need an alias is_ZnM -> is_Fp but that is it. No need for ZlM routine
because FpM_gauss use Flm_gauss internally.

> It's more like a deficiency, then. (-: Yes, I'm joking.)
> 
> > OK, I have made some test and I agree that matadjoint is very slow,
> > even compared to plain RgM_inv without the is_FpM otpimisation.
> 
> I would be interested in the numbers, if you still have them at hand.
> What I noticed is that matdet() (when using slow Gaussian elimination)
> is still usable on a slow machine, when matadjoint(,1) wouldn't work for
> me even on a really fast machine.  (It was very slow, but the worse
> thing was that it would need way too much RAM.)

Well, it is a O(n^4) algorithm instead of O(n^3) so it is not
going to be fast, but I think the problem is that RgX_RgM_eval_i is very
inefficient and waste stack space (we should use Brent&Kung).

> Somewhat unrelatedly, do you have a reference describing the Berkowitz
> algorithm and why it works?

Ask Karim...

> > There is always a cost (in storage/ memory/ testing/ maintainance/
> > consistentcy/ etc.) associate to extra code in PARI. This has to be
> > balanced with use case.
> 
> Yes, and it is not in my interest at all to send you a non-maintainable
> mess.  My understanding of your concern is that the main problem right
> now is that the amount of code may not really warrant the admittedly
> rather rare use case?

Yes, and consistency. It is always annoying to have to draft exception in
the documentation.

> > > > It would be nice if the whole strategy was abstracted to be reused in
> > > > similar
> > > > situation.
> > >
>
> On the other hand, when the modulus is detected to be non-invertible, it
> might possibly be easier to just factor the whole thing instead of doing
> the splitting trick in this case, yes.  (Although it seems to me that
> it would need probably roughly the same amount of code.)

There are other situations than matrix inversion anyway. Some operation
on ideals come to mind.
 
> > > > Do you think there is really a need for 'ZlM' ?  In that range, factoring
> > > > the modulus is trivial, and that would reduce the size of the patch
> > > > somewhat.
> > >
> > > For Flm_gauss() the problem could indeed be solved differently, but
> > > I'm not sure it's so much easier, especially since the case of p having
> > > square factors needs to be considered, and the elimination routines
> > > fixed accordingly in any case.
> > 
> > This is why I suggest to use Hensel lifting.
> 
> Ok, that can be done.
> 
> > PS:
> > Please use Fp_invsafe() (and implement a function Fl_invsafe for the Fl case).
> 
> I can do that. It means an additional duplicate gcd()-call in the
> splitting case; I don't like the look of that too much, but it can be
> done.

Feel free to use invmod() directly instead of Fp_invsafe in FpM_gauss.

> So here is my question:  What should I do to move towards resolving this
> in the best possible way?  Obviously I can get rid of
> find_coprime_factors() and fix the other issues with my last patch, and
> implement the lifting in the m = p^t case, instead of relying on the
> search for a better pivot in FpM_gaussafe().

Start by a patch that implement FpM_gausssafe, and then move on.

> But how is it integrated? I.e., how do I name the functions and where
> can I hook them in?

ZnM_inv sound fine. I will deal with the is_ZnM alias when required.

Cheers,
Bill.