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!

**Re: Flat colour in RegionPlot; millions of little triangles**

**Re: Number of interval Intersections for a large number**

**Re: Dependence of precision on execution speed of Inverse**

**RE: Dependence of precision on execution speed of Inverse**