Re: Dependence of precision on execution speed of Inverse
- To: mathgroup at smc.vnet.net
- Subject: [mg81708] Re: Dependence of precision on execution speed of Inverse
- From: "Andrew Moylan" <andrew.j.moylan at gmail.com>
- Date: Tue, 2 Oct 2007 05:30:55 -0400 (EDT)
On Sep 29, 4:28 pm, d... at wolfram.com wrote:
> On Sep 28, 1:06 am, "Andrew Moylan" <andrew.j.moy... at gmail.com> wrote:
>
>
>
> > m1=RandomReal[1,{50,50}];
> > m2=SetPrecision[m1,14];
> > Do[Inverse[m1];,{1500}]//Timing
> > Do[Inverse[m2];,{10}]//Timing
>
> > On my machine, both Timing results are about the same, indicating that
> > Inverse (of 50x50 matrices) is about 150 times faster for
machine-precision
> > numbers than for arbitrary precision numbers (of precision ~14). This
factor
> > of 150 seems large. Does Mathematica employ an Inverse algorithm that is
> > optimised for Mathematica's arbitrary-precision numbers?
>
> > Notes:
> > * Changing the precision of m2 from 14 to e.g. 17 makes little
difference.
> > * Calling Developer`FromPackedArray[] on m1 makes little difference.
> > * Calling LinearSolve on Inverse yields somewhat different results, but
> > still shows a difference in execution time of a factor of order 100.
>
> Your m1 is a packed array of machine numbers. In contrast, m2 uses
> software-implemented arbitrary precision approximate numbers. All this
> you knew; I'm just summarizing the situation.
>
> When you take Inverse[m1] you are using, under the hood, Lapack with
> optimized BLAS (Basic Linear Algebra Subroutines) code. With arbitrary
> precision, even set to a low value, there is no packing, hence no
> locality of reference for purposes of optimizing any sort of BLAS.
> That alone would force it to be much slower. Add to it the overhead of
> software implementation of the underlying arithmetic and you begin to
> see why it is so much the worse in speed.
>
> Actually if you compare with direct use of Dot you will see that
> Inverse behaves etter, relatively speaking, for the arbitrary
> precision matrix case. I'll set up a slightly different example.
>
> m1 = RandomReal[1,{100,100}];
> m2 = SetPrecision[m1,14];
> m1trans = Transpose[m1];
> m2trans = Transpose[m2];
>
> In[39]:= Timing[Do [m1.m1, {1000}];]
> Out[39]= {0.456028, Null}
>
> In[40]:= Timing[Do [m2.m2, {2}];]
> Out[40]= {2.59616, Null}
>
> So the machine arithmetic case was perhaps 3000 times faster. But we
> have internal code for handling dot products using fixed precision and
> thus bypassing issues of precision estimation speed. The test below
> shows a factor of 5 or so improvement over use of System context Dot.
>
> In[55]:= Timing[Do[
> LinearAlgebra`BLAS`GEMM["N", "N", 1, m2, m2, 0, mat], {10}]]
> Out[55]= {2.67617, Null}
>
> Beyond this, arbitrary precision Inverse is clearly gaining another
> factor of three or so relative to the machine arithmetic version,
> though offhand I know not how.
>
> All this said, there may be ways we could improve further on
> efficiency of the BLAS variants for arbitrary precision arithmetic. So
> I do not know whether the factor of 150 is higher than need be.
>
> Daniel Lichtblau
> Wolfram Research
Ah I see, thanks for your reply Daniel. So the arbitrary precision numbers
in Mathematica are each something like a struct with {a pointer to the
arbitrary precision mantissa data, the exponent, the length of the mantissa
data}; and because each element in a matrix could in principle have a
*different* arbitrary precision, there's no way to pack the array into a
contiguous lump of memory. So there's no way around de-referencing a lot of
pointers.
But Daniel, would you agree that for (hypothetical) *fixed* precision
(across the whole matrix) non-machine-precision matrices of numbers, it
*would* be possible to create the analogue of packed arrays and therefore
make optimised routines to run on them, analogous to the routines that
currently operate on packed arrays of machine-precision numbers? I guess one
can write a code (in C++ or such) that can invert non-machine-precision
matrices (and do other operations for which Mathematica employs packed
arrays only for machine precision numbers) tens of times faster than
Mathematica can, by combining something like the GNU multiple precision
library (http://gmplib.org) with BLAS-like linear algebra code.
I don't know how whether the efficiency of linear algebra at
higher-than-machine-precision affects many users, but it has come up in my
application so I am quite interested in the possibilities!