I am trying to write a function in Scheme which will calculate the rank of a given matrix (in the form of a list of lists) using Gaussian Elimination method.
I'll explain what I am trying to do in the code and would appreciate if anyone could tell me what I did wrong \ how to correct it.
In iter-rows function, I iterate over the rows and perform Gaussian Elimination on each row to make the matrix into row echelon form:
below-pivot will basically perform the elimination step below the pivot in the given row.
factor-rec here I calculate the factor needed to multiply a row during the elimination process.
elimination function performs the elimination itself over the rest of the rows below the pivot.
As for the update-matrix function, it will subtract a multiple of one row from another.
Finally, count all non-zero rows to find the rank.
Since lists are immutable I am passing their value with every call.
Here's the code for reference:
#lang racket
(define (rank matrix)
(define (iter-rows index rows-len matrix)
(if (= index rows-len)
matrix
(let ((pivot (list-ref (list-ref matrix index) index)))
(if (= pivot 0)
(iter-rows (+ index 1) rows-len matrix)
(begin
(below-pivot (+ index 1) rows-len matrix pivot index)
(iter-rows (+ index 1) rows-len matrix))))))
(define (below-pivot from-index rows-len matrix pivot i)
(if (= from-index rows-len)
matrix
(begin
(factor-rec (+ from-index 1) rows-len matrix i pivot)
(below-pivot (+ from-index 1) rows-len matrix pivot i))))
(define (factor-rec j rows-len matrix i pivot)
(if (= j rows-len)
matrix
(begin
(let ([factor (/ (list-ref (list-ref matrix j) i) pivot)])
(elimination i i rows-len matrix j factor)))))
(define (elimination k i rows-len matrix j factor)
(if (= k rows-len)
matrix
(begin
(update-matrix matrix j k i factor)
(elimination (+ k 1) i rows-len matrix j factor))))
(define (update-matrix matrix j k i factor)
(let ((new-cdr (- (list-ref (list-ref matrix j) k)
(* factor (list-ref (list-ref matrix i) k)))))
(list (list-ref matrix i)
(cons new-cdr (list-tail (list-ref matrix j) k)))))
(let ((rows-len (length matrix)))
(iter-rows 0 rows-len matrix))
(define (count-nonzero-rows matrix)
(let loop ((matrix matrix) (rank- 0))
(cond
((null? matrix) rank-)
((not (all-zeroes? (car matrix)))
(loop (cdr matrix) (+ rank- 1)))
(else (loop (cdr matrix) rank-)))))
(let* ((rank (count-nonzero-rows matrix)))
rank))
(define (all-zeroes? row)
(for/and ([x row]) (= x 0)))
I'm new to Scheme so would appreciate any other tips as well.
Thanks.
I think you want something like this:
#lang racket
(define (extract matrix i j)
(vector-ref (vector-ref matrix i) j))
(define (swap-rows! matrix i j)
(let ((tmp (vector-ref matrix i)))
(vector-set! matrix i (vector-ref matrix j))
(vector-set! matrix j tmp)))
(define (scale-row! matrix i factor)
(define (vector-scale vec factor)
(vector-map (curry * factor) vec))
(let* ((row (vector-ref matrix i))
(scaled-row (vector-scale row factor)))
(vector-set! matrix i scaled-row)))
(define (eliminate-rows! matrix i j pivot)
(define (eliminate-row! matrix i j pivot)
(when (< i (vector-length matrix))
(let* ((element-i-j (extract matrix i j))
(factor (/ element-i-j pivot)))
(scale-row! matrix i factor)
(eliminate-row! matrix (add1 i) j pivot))))
(eliminate-row! matrix (add1 i) j pivot))
(define (gaussian-rank matrix (i 0) (j 0) (rank 0))
(cond ((>= i (vector-length matrix)) rank)
((zero? (extract matrix i j)) (gaussian-rank matrix (add1 i) j rank))
(else
(swap-rows! matrix i rank)
(let ((pivot (extract matrix rank j)))
(when (not (zero? pivot))
(scale-row! matrix rank (/ 1 pivot))
(eliminate-rows! matrix rank j pivot))
(gaussian-rank matrix (add1 i) (add1 j) (add1 rank))))))
(define mat (vector #(1 2 3) #(0 1 2) #(0 0 1)))
(gaussian-rank mat) ;;=> 1
The "matrix" should be a vector. (Probably one should introduce a struct).
extract returns the element in a matrix in the i-th row and j-th column.
swap-rows! takes a matrix and rows the i-th and j-th row.
scale-row! scales a row multiplying each of the i-th row element by a factor.
eliminate-rows! uses these two functions to eliminat-through a matrix.
And these operations are used to calculat the gaussian rank.
Perhaps more efficient is this version:
(define (get-row matrix i)
(vector-ref matrix i))
(define (get-element matrix i j)
(vector-ref (vector-ref matrix i) j))
(define (swap-rows! matrix i j)
(let ((tmp (vector-ref matrix i)))
(vector-set! matrix i (vector-ref matrix j))
(vector-set! matrix j tmp)))
(define-syntax while
(syntax-rules () ((_ pred b1 ...) (let loop () (when pred b1 ... (loop))))))
(define (gaussian-rank matrix)
(let* ((nrows (vector-length matrix))
(ncols (if (zero? nrows) 0 (vector-length (vector-ref matrix 0))))
(rank 0))
(for ([ncol (in-range ncols)])
(let ([pivot-row rank])
(while (and (< pivot-row nrows)
(zero? (get-element matrix pivot-row ncol)))
(set! pivot-row (add1 pivot-row)))
(unless (= pivot-row nrows)
(swap-rows! matrix pivot-row rank)
(let ([pivot-value (get-element matrix rank ncol)])
(for ([nrow (in-range (add1 rank) nrows)])
(let* ([factor (/ (get-element matrix nrow ncol) pivot-value)]
[scaled-diffed-row (vector-map (lambda (x y) (- x (* factor y)))
(get-row matrix nrow)
(get-row matrix rank))])
(vector-set! matrix nrow scaled-diffed-row))))
(set! rank (add1 rank)))))
rank))
Here, the same processing is done but working preferably with loops and thus perhaps faster.
Call it by:
> (define mat (vector #(1 2 3) #(4 5 6) #(7 8 9)))
> (gaussian-rank mat)
2
> mat
'#(#(1 2 3) #(0 -3 -6) #(0 0 0))
the syntax-rule takes:
(while (and (< pivot-row nrows)
(zero? (get-element matrix pivot-row ncol)))
(set! pivot-row (add1 pivot-row)))
and matches
`_`: while
`pred`: (and (< pivot-row nrows)
(zero? (get-element matrix pivot-row ncol)))
`b1`: (set! pivot-row (add1 pivot-row))
`...`: in this case nothing
and transforms this - while nothing is yet evaluated - to:
(let loop ()
(when (and (< pivot-row nrows)
(zero? (get-element matrix pivot-row ncol)))
(set! pivot-row (add1 pivot-row))
(loop))
and then evalutes this generated code in-situ (in exactly the environment/place where the while expression was invoked).
This code transformation is done in compile time.
This is how a lisp macro works.