[racket] math/matrix
That's great!
The question is now how to automate this sort of thing.
/Jens Axel
2014-05-13 1:39 GMT+02:00 Neil Toronto <neil.toronto at gmail.com>:
> When I change it to operate on (Vectorof FlVector) instead of (Vectorof
> (Vectorof Flonum)), I get this:
>
> cpu time: 996 real time: 995 gc time: 22
> 1.0000000000009335
> cpu time: 15387 real time: 15384 gc time: 13006
> 1.0000000000009335
> cpu time: 1057 real time: 1056 gc time: 85
> 1.0000000000009335
> cpu time: 11514 real time: 11510 gc time: 9097
> 1.0000000000009335
> cpu time: 1079 real time: 1079 gc time: 100
> 1.0000000000009335
> cpu time: 15425 real time: 15426 gc time: 13072
> 1.0000000000009335
>
> When I use `racket/unsafe/ops` instead of `math/private/unsafe`, I get this:
>
> cpu time: 591 real time: 591 gc time: 17
> 1.0000000000016622
> cpu time: 11514 real time: 11509 gc time: 9195
> 1.0000000000016622
> cpu time: 604 real time: 604 gc time: 31
> 1.0000000000016622
> cpu time: 15739 real time: 15737 gc time: 13358
> 1.0000000000016622
> cpu time: 596 real time: 595 gc time: 24
> 1.0000000000016622
> cpu time: 11498 real time: 11493 gc time: 9154
> 1.0000000000016622
>
> Racket's floating-point math is fast if the flonums aren't allocated
> separately on the heap.
>
> Neil ⊥
>
>
> On 05/12/2014 02:17 PM, Jens Axel Søgaard wrote:
>>
>> Hi Eric,
>>
>> You were absolute right. The version below cuts the time in half.
>> It is mostly cut and paste from existing functions and removing
>> non-Flonum cases.
>>
>> /Jens Axel
>>
>> #lang typed/racket
>> (require math/matrix
>> math/array
>> math/private/matrix/utils
>> math/private/vector/vector-mutate
>> math/private/unsafe
>> (only-in racket/unsafe/ops unsafe-fl/)
>> racket/fixnum
>> racket/flonum
>> racket/list)
>>
>> (define-type Pivoting (U 'first 'partial))
>>
>> (: flonum-matrix-gauss-elim
>> (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
>> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof
>> Index)))
>> ((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof
>> Index)))
>> ((Matrix Flonum) Any Any Pivoting -> (Values (Matrix
>> Flonum) (Listof Index)))))
>> (define (flonum-matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]
>> [pivoting 'partial])
>> (define-values (m n) (matrix-shape M))
>> (define rows (matrix->vector* M))
>> (let loop ([#{i : Nonnegative-Fixnum} 0]
>> [#{j : Nonnegative-Fixnum} 0]
>> [#{without-pivot : (Listof Index)} empty])
>> (cond
>> [(j . fx>= . n)
>> (values (vector*->matrix rows)
>> (reverse without-pivot))]
>> [(i . fx>= . m)
>> (values (vector*->matrix rows)
>> ;; None of the rest of the columns can have pivots
>> (let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot
>> without-pivot])
>> (cond [(j . fx< . n) (loop (fx+ j 1) (cons j
>> without-pivot))]
>> [else (reverse without-pivot)])))]
>> [else
>> (define-values (p pivot)
>> (case pivoting
>> [(partial) (find-partial-pivot rows m i j)]
>> [(first) (find-first-pivot rows m i j)]))
>> (cond
>> [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
>> [else
>> ;; Swap pivot row with current
>> (vector-swap! rows i p)
>> ;; Possibly unitize the new current row
>> (let ([pivot (if unitize-pivot?
>> (begin (vector-scale! (unsafe-vector-ref rows
>> i)
>> (unsafe-fl/ 1. pivot))
>> (unsafe-fl/ pivot pivot))
>> pivot)])
>> (flonum-elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
>> (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
>>
>> (: flonum-elim-rows!
>> ((Vectorof (Vectorof Flonum)) Index Index Index Flonum
>> Nonnegative-Fixnum -> Void))
>> (define (flonum-elim-rows! rows m i j pivot start)
>> (define row_i (unsafe-vector-ref rows i))
>> (let loop ([#{l : Nonnegative-Fixnum} start])
>> (when (l . fx< . m)
>> (unless (l . fx= . i)
>> (define row_l (unsafe-vector-ref rows l))
>> (define x_lj (unsafe-vector-ref row_l j))
>> (unless (= x_lj 0)
>> (flonum-vector-scaled-add! row_l row_i (fl* -1. (fl/ x_lj
>> pivot)) j)
>> ;; Make sure the element below the pivot is zero
>> (unsafe-vector-set! row_l j (- x_lj x_lj))))
>> (loop (fx+ l 1)))))
>>
>>
>> (: flonum-matrix-solve
>> (All (A) (case->
>> ((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
>> ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix
>> Flonum))))))
>> (define flonum-matrix-solve
>> (case-lambda
>> [(M B) (flonum-matrix-solve
>> M B (λ () (raise-argument-error 'flonum-matrix-solve
>> "matrix-invertible?" 0 M B)))]
>> [(M B fail)
>> (define m (square-matrix-size M))
>> (define-values (s t) (matrix-shape B))
>> (cond [(= m s)
>> (define-values (IX wps)
>> (parameterize ([array-strictness #f])
>> (flonum-matrix-gauss-elim (matrix-augment (list M B)) #t
>> #t)))
>> (cond [(and (not (empty? wps)) (= (first wps) m))
>> (submatrix IX (::) (:: m #f))]
>> [else (fail)])]
>> [else
>> (error 'flonum-matrix-solve
>> "matrices must have the same number of rows; given ~e
>> and ~e"
>> M B)])]))
>>
>> (define-syntax-rule (flonum-vector-generic-scaled-add! vs0-expr
>> vs1-expr v-expr start-expr + *)
>> (let* ([vs0 vs0-expr]
>> [vs1 vs1-expr]
>> [v v-expr]
>> [n (fxmin (vector-length vs0) (vector-length vs1))])
>> (let loop ([#{i : Nonnegative-Fixnum} (fxmin start-expr n)])
>> (if (i . fx< . n)
>> (begin (unsafe-vector-set! vs0 i (+ (unsafe-vector-ref vs0 i)
>> (* (unsafe-vector-ref vs1
>> i) v)))
>> (loop (fx+ i 1)))
>> (void)))))
>>
>> (: flonum-vector-scaled-add!
>> (case-> ((Vectorof Flonum) (Vectorof Flonum) Flonum -> Void)
>> ((Vectorof Flonum) (Vectorof Flonum) Flonum Index -> Void)))
>> (define (flonum-vector-scaled-add! vs0 vs1 s [start 0])
>> (flonum-vector-generic-scaled-add! vs0 vs1 s start + *))
>>
>> (: mx Index)
>> (define mx 600)
>>
>> (: r (Index Index -> Flonum))
>> (define (r i j) (random))
>>
>> (: A : (Matrix Flonum))
>> (define A (build-matrix mx mx r))
>>
>> (: sum : Integer Integer -> Flonum)
>> (define (sum i n)
>> (let loop ((j 0) (acc 0.0))
>> (if (>= j mx) acc
>> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>>
>> (: b : (Matrix Flonum))
>> (define b (build-matrix mx 1 sum))
>>
>> (time
>> (let [(m (flonum-matrix-solve A b))]
>> (matrix-ref m 0 0)))
>> (time
>> (let [(m (matrix-solve A b))]
>> (matrix-ref m 0 0)))
>>
>> (time
>> (let [(m (flonum-matrix-solve A b))]
>> (matrix-ref m 0 0)))
>> (time
>> (let [(m (matrix-solve A b))]
>> (matrix-ref m 0 0)))
>>
>> (time
>> (let [(m (flonum-matrix-solve A b))]
>> (matrix-ref m 0 0)))
>> (time
>> (let [(m (matrix-solve A b))]
>> (matrix-ref m 0 0)))
>>
>> /Jens Axel
>>
>>
>> 2014-05-11 23:26 GMT+02:00 Eric Dobson <eric.n.dobson at gmail.com>:
>>>
>>> Where is the time spent in the algorithm? I assume that most of it is
>>> in the matrix manipulation work not the orchestration of finding a
>>> pivot and reducing based on that. I.e. `elim-rows!` is the expensive
>>> part. Given that you only specialized the outer part of the loop, I
>>> wouldn't expect huge performance changes.
>>>
>>> On Sun, May 11, 2014 at 2:13 PM, Jens Axel Søgaard
>>> <jensaxel at soegaard.net> wrote:
>>>>
>>>> I tried restricting the matrix-solve and matrix-gauss-elim to (Matrix
>>>> Flonum).
>>>> I can't observe a change in the timings.
>>>>
>>>> #lang typed/racket
>>>> (require math/matrix
>>>> math/array
>>>> math/private/matrix/utils
>>>> math/private/vector/vector-mutate
>>>> math/private/unsafe
>>>> (only-in racket/unsafe/ops unsafe-fl/)
>>>> racket/fixnum
>>>> racket/list)
>>>>
>>>> (define-type Pivoting (U 'first 'partial))
>>>>
>>>> (: flonum-matrix-gauss-elim
>>>> (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
>>>> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof
>>>> Index)))
>>>> ((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof
>>>> Index)))
>>>> ((Matrix Flonum) Any Any Pivoting -> (Values (Matrix
>>>> Flonum) (Listof Index)))))
>>>> (define (flonum-matrix-gauss-elim M [jordan? #f] [unitize-pivot? #f]
>>>> [pivoting 'partial])
>>>> (define-values (m n) (matrix-shape M))
>>>> (define rows (matrix->vector* M))
>>>> (let loop ([#{i : Nonnegative-Fixnum} 0]
>>>> [#{j : Nonnegative-Fixnum} 0]
>>>> [#{without-pivot : (Listof Index)} empty])
>>>> (cond
>>>> [(j . fx>= . n)
>>>> (values (vector*->matrix rows)
>>>> (reverse without-pivot))]
>>>> [(i . fx>= . m)
>>>> (values (vector*->matrix rows)
>>>> ;; None of the rest of the columns can have pivots
>>>> (let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot
>>>> without-pivot])
>>>> (cond [(j . fx< . n) (loop (fx+ j 1) (cons j
>>>> without-pivot))]
>>>> [else (reverse without-pivot)])))]
>>>> [else
>>>> (define-values (p pivot)
>>>> (case pivoting
>>>> [(partial) (find-partial-pivot rows m i j)]
>>>> [(first) (find-first-pivot rows m i j)]))
>>>> (cond
>>>> [(zero? pivot) (loop i (fx+ j 1) (cons j without-pivot))]
>>>> [else
>>>> ;; Swap pivot row with current
>>>> (vector-swap! rows i p)
>>>> ;; Possibly unitize the new current row
>>>> (let ([pivot (if unitize-pivot?
>>>> (begin (vector-scale! (unsafe-vector-ref
>>>> rows i)
>>>> (unsafe-fl/ 1.
>>>> pivot))
>>>> (unsafe-fl/ pivot pivot))
>>>> pivot)])
>>>> (elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
>>>> (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
>>>>
>>>> (: flonum-matrix-solve
>>>> (All (A) (case->
>>>> ((Matrix Flonum) (Matrix Flonum) -> (Matrix
>>>> Flonum))
>>>> ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix
>>>> Flonum))))))
>>>> (define flonum-matrix-solve
>>>> (case-lambda
>>>> [(M B) (flonum-matrix-solve
>>>> M B (λ () (raise-argument-error 'flonum-matrix-solve
>>>> "matrix-invertible?" 0 M B)))]
>>>> [(M B fail)
>>>> (define m (square-matrix-size M))
>>>> (define-values (s t) (matrix-shape B))
>>>> (cond [(= m s)
>>>> (define-values (IX wps)
>>>> (parameterize ([array-strictness #f])
>>>> (flonum-matrix-gauss-elim (matrix-augment (list M B))
>>>> #t #t)))
>>>> (cond [(and (not (empty? wps)) (= (first wps) m))
>>>> (submatrix IX (::) (:: m #f))]
>>>> [else (fail)])]
>>>> [else
>>>> (error 'flonum-matrix-solve
>>>> "matrices must have the same number of rows; given
>>>> ~e and ~e"
>>>> M B)])]))
>>>>
>>>> (: mx Index)
>>>> (define mx 600)
>>>>
>>>> (: r (Index Index -> Flonum))
>>>> (define (r i j) (random))
>>>>
>>>> (: A : (Matrix Flonum))
>>>> (define A (build-matrix mx mx r))
>>>>
>>>> (: sum : Integer Integer -> Flonum)
>>>> (define (sum i n)
>>>> (let loop ((j 0) (acc 0.0))
>>>> (if (>= j mx) acc
>>>> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>>>>
>>>> (: b : (Matrix Flonum))
>>>> (define b (build-matrix mx 1 sum))
>>>>
>>>> (time
>>>> (let [(m (flonum-matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>> (time
>>>> (let [(m (matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>>
>>>> (time
>>>> (let [(m (flonum-matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>> (time
>>>> (let [(m (matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>>
>>>> (time
>>>> (let [(m (flonum-matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>> (time
>>>> (let [(m (matrix-solve A b))]
>>>> (matrix-ref m 0 0)))
>>>>
>>>> 2014-05-11 21:48 GMT+02:00 Neil Toronto <neil.toronto at gmail.com>:
>>>>>
>>>>> The garbage collection time is probably from cleaning up boxed flonums,
>>>>> and
>>>>> possibly intermediate vectors. If so, a separate implementation of
>>>>> Gaussian
>>>>> elimination for the FlArray type would cut the GC time to nearly zero.
>>>>>
>>>>> Neil ⊥
>>>>>
>>>>>
>>>>> On 05/11/2014 01:36 PM, Jens Axel Søgaard wrote:
>>>>>>
>>>>>>
>>>>>> Or ... you could take a look at
>>>>>>
>>>>>>
>>>>>>
>>>>>> https://github.com/plt/racket/blob/master/pkgs/math-pkgs/math-lib/math/private/matrix/matrix-gauss-elim.rkt
>>>>>>
>>>>>> at see if something can be improved.
>>>>>>
>>>>>> /Jens Axel
>>>>>>
>>>>>>
>>>>>> 2014-05-11 21:30 GMT+02:00 Jens Axel Søgaard <jensaxel at soegaard.net>:
>>>>>>>
>>>>>>>
>>>>>>> Hi Eduardo,
>>>>>>>
>>>>>>> The math/matrix library uses the arrays from math/array to represent
>>>>>>> matrices.
>>>>>>>
>>>>>>> If you want to try the same representation as Bigloo, you could try
>>>>>>> Will
>>>>>>> Farr's
>>>>>>> matrix library:
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> http://planet.racket-lang.org/package-source/wmfarr/simple-matrix.plt/1/1/planet-docs/simple-matrix/index.html
>>>>>>>
>>>>>>> I am interested in hearing the results.
>>>>>>>
>>>>>>> /Jens Axel
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> 2014-05-11 21:18 GMT+02:00 Eduardo Costa <edu500ac at gmail.com>:
>>>>>>>>
>>>>>>>>
>>>>>>>> What is bothering me is the time Racket is spending in garbage
>>>>>>>> collection.
>>>>>>>>
>>>>>>>> ~/wrk/scm/rkt/matrix$ racket matrix.rkt
>>>>>>>> 0.9999999999967226
>>>>>>>> cpu time: 61416 real time: 61214 gc time: 32164
>>>>>>>>
>>>>>>>> If I am reading the output correctly, Racket is spending 32 seconds
>>>>>>>> out
>>>>>>>> of
>>>>>>>> 61 seconds in garbage collection.
>>>>>>>>
>>>>>>>> I am following Junia Magellan's computer language comparison and
>>>>>>>> I
>>>>>>>> cannot
>>>>>>>> understand why Racket needs the garbage collector for doing Gaussian
>>>>>>>> elimination. In a slow Compaq/HP machine, solving a system of 800
>>>>>>>> linear
>>>>>>>> equations takes 17.3 seconds in Bigloo, but requires 58 seconds in
>>>>>>>> Racket,
>>>>>>>> even after removing the building of the linear system from
>>>>>>>> consideration.
>>>>>>>> Common Lisp is also much faster than Racket in processing arrays. I
>>>>>>>> would
>>>>>>>> like to point out that Racket is very fast in general. The only
>>>>>>>> occasion
>>>>>>>> that it lags badly behind Common Lisp and Bigloo is when one needs
>>>>>>>> to
>>>>>>>> deal
>>>>>>>> with arrays.
>>>>>>>>
>>>>>>>> Basically, Junia is using Rasch method to measure certain latent
>>>>>>>> traits
>>>>>>>> of
>>>>>>>> computer languages, like productivity and coaching time. In any
>>>>>>>> case,
>>>>>>>> she
>>>>>>>> needs to do a lot of matrix calculations to invert the Rasch model.
>>>>>>>> Since
>>>>>>>> Bigloo works with homogeneous vectors, she wrote a few macros to
>>>>>>>> access
>>>>>>>> the
>>>>>>>> elements of a matrix:
>>>>>>>>
>>>>>>>> (define (mkv n) (make-f64vector n))
>>>>>>>> (define $ f64vector-ref)
>>>>>>>> (define $! f64vector-set!)
>>>>>>>> (define len f64vector-length)
>>>>>>>>
>>>>>>>> (define-syntax $$
>>>>>>>> (syntax-rules ()
>>>>>>>> (($$ m i j) (f64vector-ref (vector-ref m i) j))))
>>>>>>>>
>>>>>>>> (define-syntax $$!
>>>>>>>> (syntax-rules ()
>>>>>>>> (($$! matrix row column value)
>>>>>>>> ($! (vector-ref matrix row) column value))))
>>>>>>>>
>>>>>>>> I wonder whether homogeneous vectors would speed up Racket. In the
>>>>>>>> same
>>>>>>>> computer that Racket takes 80 seconds to build and invert a system
>>>>>>>> of
>>>>>>>> equations, Bigloo takes 17.3 seconds, as I told before. Common Lisp
>>>>>>>> is
>>>>>>>> even
>>>>>>>> faster. However, if one subtracts the gc time from Racket's total
>>>>>>>> time,
>>>>>>>> the
>>>>>>>> result comes quite close to Common Lisp or Bigloo.
>>>>>>>>
>>>>>>>> ~/wrk/bgl$ bigloo -Obench bigmat.scm -o big
>>>>>>>> ~/wrk/bgl$ time ./big
>>>>>>>> 0.9999999999965746 1.000000000000774 0.9999999999993039
>>>>>>>> 0.9999999999982576
>>>>>>>> 1.000000000007648 0.999999999996588
>>>>>>>>
>>>>>>>> real 0m17.423s
>>>>>>>> user 0m17.384s
>>>>>>>> sys 0m0.032s
>>>>>>>> ~/wrk/bgl$
>>>>>>>>
>>>>>>>> Well, bigloo may perform global optimizations, but Common Lisp
>>>>>>>> doesn't.
>>>>>>>> When
>>>>>>>> one is not dealing with matrices, Racket is faster than Common Lisp.
>>>>>>>> I
>>>>>>>> hope
>>>>>>>> you can tell me how to rewrite the program in order to avoid garbage
>>>>>>>> collection.
>>>>>>>>
>>>>>>>> By the way, you may want to know why not use Bigloo or Common Lisp
>>>>>>>> to
>>>>>>>> invert
>>>>>>>> the Rasch model. The problem is that Junia and her co-workers are
>>>>>>>> using
>>>>>>>> hosting services that do not give access to the server or to the
>>>>>>>> jailshell.
>>>>>>>> Since Bigloo requires gcc based compilation, Junia discarded it
>>>>>>>> right
>>>>>>>> away.
>>>>>>>> Not long ago, the hosting service stopped responding to the sbcl
>>>>>>>> Common
>>>>>>>> Lisp
>>>>>>>> compiler for reasons that I cannot fathom. Although Racket 6.0
>>>>>>>> stopped
>>>>>>>> working too, Racket 6.0.1 is working fine. This left Junia, her
>>>>>>>> co-workers
>>>>>>>> and students with Racket as their sole option. As for myself, I am
>>>>>>>> just
>>>>>>>> curious.
>>>>>>>>
>>>>>>>>
>>>>>>>> 2014-05-11 6:23 GMT-03:00 Jens Axel Søgaard <jensaxel at soegaard.net>:
>>>>>>>>
>>>>>>>>> 2014-05-11 6:09 GMT+02:00 Eduardo Costa <edu500ac at gmail.com>:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> The documentation says that one should expect typed/racket to be
>>>>>>>>>> faster
>>>>>>>>>> than
>>>>>>>>>> racket. I tested the math/matrix library and it seems to be almost
>>>>>>>>>> as
>>>>>>>>>> slow
>>>>>>>>>> in typed/racket as in racket.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> What was (is?) slow was a call in an untyped module A to a function
>>>>>>>>> exported
>>>>>>>>> from a typed module B. The functions in B must check at runtime
>>>>>>>>> that
>>>>>>>>> the values coming from A are of the correct type. If the A was
>>>>>>>>> written
>>>>>>>>> in Typed Racket, the types would be known at compile time.
>>>>>>>>>
>>>>>>>>> Here math/matrix is written in Typed Racket, so if you are writing
>>>>>>>>> an
>>>>>>>>> untyped module, you will in general want to minimize the use
>>>>>>>>> of,say,
>>>>>>>>> maxtrix-ref. Instead operations that works on entire matrices or
>>>>>>>>> row/columns are preferred.
>>>>>>>>>
>>>>>>>>>> (: sum : Integer Integer -> Flonum)
>>>>>>>>>> (define (sum i n)
>>>>>>>>>> (let loop ((j 0) (acc 0.0))
>>>>>>>>>> (if (>= j mx) acc
>>>>>>>>>> (loop (+ j 1) (+ acc (matrix-ref A i j))) )))
>>>>>>>>>>
>>>>>>>>>> (: b : (Matrix Flonum))
>>>>>>>>>> (define b (build-matrix mx 1 sum))
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> The matrix b contains the sums of each row in the matrix.
>>>>>>>>> Since matrices are a subset of arrays, you can use array-axis-sum,
>>>>>>>>> which computes sum along a given axis (i.e. a row or a column when
>>>>>>>>> speaking of matrices).
>>>>>>>>>
>>>>>>>>> (define A (matrix [[0. 1. 2.]
>>>>>>>>> [3. 4. 5.]
>>>>>>>>> [6. 7. 8.]]))
>>>>>>>>>
>>>>>>>>>> (array-axis-sum A 1)
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> - : (Array Flonum)
>>>>>>>>> (array #[3.0 12.0 21.0])
>>>>>>>>>
>>>>>>>>> However as Eric points out, matrix-solve is an O(n^3) algorithm,
>>>>>>>>> so the majority of the time is spent in matrix-solve.
>>>>>>>>>
>>>>>>>>> Apart from finding a way to exploit the relationship between your
>>>>>>>>> matrix A and the column vector b, I see no obvious way of
>>>>>>>>> speeding up the code.
>>>>>>>>>
>>>>>>>>> Note that when you benchmark with
>>>>>>>>>
>>>>>>>>> time racket matrix.rkt
>>>>>>>>>
>>>>>>>>> you will include startup and compilation time.
>>>>>>>>> Therefore if you want to time the matrix code,
>>>>>>>>> insert a literal (time ...) call.
>>>>>>>>>
>>>>>>>>> --
>>>>>>>>> Jens Axel Søgaard
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> --
>>>>>>> Jens Axel Søgaard
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> ____________________
>>>>> Racket Users list:
>>>>> http://lists.racket-lang.org/users
>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> --
>>>> Jens Axel Søgaard
>>>>
>>>> ____________________
>>>> Racket Users list:
>>>> http://lists.racket-lang.org/users
>>
>>
>>
>>
>
--
--
Jens Axel Søgaard
![[logo]](/logo.png){kind=logo}