1000 Equations With 1000 Unknowns

[Dick Ricker]

FreeMat takes a very long time to solve large sets of equations using rref(). On my computer (3 GHz Pentium 4, 1 GB RAM, Windows XP, SP3, FreeMat 4.2 for Windows XP released August 15, 2013), rref() with a [1000 x 1001] matrix takes almost an hour to run. Probably because rref() is not JIT Compile friendly.

A much faster alternative is to use the multiply by inverse approach:

    AX = B                        :  equation in vector form

    (A^-1)AX = (A^-1)B       :  left multiply both sides by (A^-1)

    IX = (A^-1)B                 :  (A^-1)A = identity matrix

    X = (A^-1)B                  : IX = X

Test Program:

printf('***** [100 x 100] + [100 x 1]\n\n')

A = fix(100 * rand(100)) ;

B = fix(100 * rand(100, 1)) ;

AugA = [ A B ] ;

time1 = clock ;

Ainv = A^-1 ;

Xinv = (Ainv * B) ;

time2 = clock ;

exe_time21 = etime(time2, time1) ;

printf('Invert A, Multiply B = %5.3f Seconds\n\n', exe_time21)

time3 = clock ;

RredA = rref(AugA) ;

Xred = RredA( :, 101 ) ;

time4 = clock ;

exe_time34 = etime(time4, time3) ;

printf('Row reduce AugA, Extract X = %5.3f Seconds\n\n', exe_time34)

Xdif = Xred - Xinv ;

Xabs = abs(Xdif) ;

Xmax = max(Xabs) ;

printf('Maximum absolute difference of X vectors = %e\n\n', Xmax)

Test Results:

[100 x 100] + [100 x 1]

    invert A, multiply B : 0.016 seconds

    row reduce AugA, Extract X : 3.66 seconds

    maximum difference : 1.15e-14

[200 x 200] + [200 x 1]

    invert A, multiply B : 0.110 seconds

    row reduce AugA, Extract X : 27.3 seconds

    maximum difference : 1.12e-13

[300 x 300] + [300 x 1]

    invert A, multiply B : 0.344 seconds

    row reduce AugA, Extract X : 90.7 seconds

    maximum difference : 1.29e-13

[500 x 500] + [500 x 1]

    invert A, multiply B : 1.56 seconds

    row reduce AugA, Extract X : 414 seconds

    maximum difference : 3.19e-13

[1000 x 1000] + [1000 x 1]

    invert A, multiply B : 12.6 seconds

    row reduce AugA, Extract X : 3273 seconds

    maximum difference : 3.68e-11

The only downside of the multiply by inverse approach is that the A matrix must be invertible. That means no rows or columns of all zeros, and no rows or columns which are the same. This is not a serious restriction, because in those cases there would not be a unique solution set anyway.

I have attached my complete test program for those who would like to give it a try.

Dick Ricker

[Timothy Cyders]

Dick,

I'm curious as to why you wouldn't use Gaussian elimination with A\B to solve the system instead. My 3 GHz processor takes only several milliseconds to perform the same operation to the same level of accuracy.

Cheers,

TJ

--
You received this message because you are subscribed to the Google Groups "freemat" group.
To unsubscribe from this group and stop receiving emails from it, send an email to freemat+u...@googlegroups.com.
To post to this group, send email to fre...@googlegroups.com.
Visit this group at http://groups.google.com/group/freemat.
For more options, visit https://groups.google.com/groups/opt_out.

[Timothy Cyders]

As an addendum, what you're seeing is an essentially linear increase in the number of seconds required for the operation with the number of elements in the array passed to rref(). The single line in the program responsible for most of the time spent is line 76, where a large matrix multiplication operation takes place. You can detect these time sinks using the profiler command. You can turn it on using 'profiler on', run a script or function, and then type 'profiler list functionname' to see a breakdown of what parts of the program are consuming the most processor time (sans quotes, of course). Also, for a truncated version of your clock commands, you can use the tic and toc commands.

Cheers,

TJ

[Dick Ricker]

Thank You, TJ!

I will try A\B and also try out the profiler.

Dick Ricker

[Dick Ricker]

Hi TJ,

Here is what I found:

Left Divide (A\B) is significantly faster than (A^-1)B.  For a [1000 x 1000] + [1000 x 1] set, I get about 4.5 seconds using A\B, where (A^-1)B was taking about 12.5 seconds, and rref() was taking about 3275 seconds.  And yes, error is in the same order of magnitude for all three: e-11.

tic and toc are a nice reduction of code.

profiler shows that 91% of my program is spent on A\B

Thanks Again!

Dick Ricker

PS : I will post my test results when I get a little more time.

[Dick Ricker]

Left Divide (A\B) Test Results :

[1000 x 1000] + [1000 x 1]

     A\B : 4.34 seconds

     maximum error : 7.25e-12

[2000 x 2000] + [2000 x 1]

     A\B : 33.9 seconds

     maximum error : 6.80e-11

[3000 x 3000] + [3000 x 1]

     A\B : 112 seconds

     maximum error : 1.73e-11

[4000 x 4000] + [4000 x 1]

     A\B : 264 seconds

     maximum error : 3.91e-11

[4500 x 4500] + [4500 x 1]

     A\B : 381 seconds

     maximum error : 9.33e-10

[5000 x 5000] + [5000 x 1]

     Error: Cannot allocate enough memory to store an

     array of size 25000000

I have attached my test program for those who are interested.

Dick Ricker