[racket] math/matrix

From: Matthias Felleisen (matthias at ccs.neu.edu)
Date: Tue May 13 17:24:27 EDT 2014

On May 13, 2014, at 4:19 PM, Neil Toronto <neil.toronto at gmail.com> wrote:

> We need a predicate like
> 
>  (: flonum-matrix? (All (A) (-> (Matrix A) Boolean : (Matrix Flonum))))


I think in our world of types we could even have 

 (: flonum-matrix? (All (A) (-> (Matrix A) Boolean : (TriangularMatrix A))))

and such and then dispatch to even more special solvers. It's kind of like a number hierarchy generalization. Just a thought. 


> 
> Then `matrix-solve` could dispatch to `flmatrix-solve` and still be well-typed. We could/should do something similar for every operation for which checking flonum-ness is cheap compared to computing the result, which at least includes everything O(n^3).
> 
> One thing we should really do is get your LAPACK FFI into the math library and have `flmatrix-solve` use that, but fail over to Racket code systems that don't have LAPACK. If I remember right, it would have to transpose the data because LAPACK is column-major.
> 
> Some thoughts, in no particular order:
> 
> 1. Because of transposition and FFI overhead, there's a matrix size threshold under which we ideally should use the code below, even on systems with LAPACK installed.
> 
> 2. Because of small differences in how it finds pivots, LAPACK's solver can return slightly different results. Should we worry about that at all?
> 
> 3. A design decision: if a matrix contains just one flonum, should we convert it to (Matrix Flonum) and solve it quickly with `flmatrix-solve`, or use the current `matrix-solve` to preserve some of its exactness?
> 
> I lean toward regarding a matrix with one flonum as a flonum matrix. It's definitely easier to write library code for, and would make it easier for users to predict when a result entry will be exact. Currently, we have this somewhat confusing situation, in which how pivots are chosen determines which result entries are exact:
> 
> > (matrix-row-echelon (matrix ([1 2 3] [4.0 5 4])) #t #t 'first)
> (mutable-array #[#[1 0 -2.333333333333333]
>                 #[-0.0 1.0 2.6666666666666665]])
> 
> > (matrix-row-echelon (matrix ([1 2 3] [4.0 5 4])) #t #t 'partial)
> (mutable-array #[#[1.0 0.0 -2.333333333333333]
>                 #[0 1.0 2.6666666666666665]])
> 
> I doubt this has caused problems for anyone, but it bothers me a little.
> 
> Neil ⊥
> 
> On 05/13/2014 12:45 PM, Jens Axel Søgaard wrote:
>> 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
>>>> 
>>>> 
>>>> 
>>>> 
>>> 
>> 
>> 
>> 
> 
> ____________________
> Racket Users list:
> http://lists.racket-lang.org/users



Posted on the users mailing list.