There is a special place in hell for all of those library writers who
implemented rank as the number of indices needed to access
elements of a multi-dimensional array (that’s ‘tensor’ for you!). Even
common lisp does it! What’s up with that?!
So, today I am going to implement the correct rank in common lisp.
The implementation below is based on the Gram-Schmidt
process. But first, I will need few things from gsll
(ql:quickload :gsll)(GSLL)and here is the implementation:
(defun linalg-rank (xs &optional (n 0) (epsilon 5d-3))
(destructuring-bind (a b) (array-dimensions xs)
(if (> a b)
(linalg-rank (grid:transpose xs) 0)
(let* ((ys (make-array b :displaced-to xs :displaced-index-offset 0))
(zs (grid:normalize ys))
(d (grid:norm ys))
(m (if (< d epsilon) n (1+ n)))
(vs (make-array (list (1- a) b) :displaced-to xs :displaced-index-offset b)))
(cond
((= a 1) m)
((< d epsilon) (linalg-rank vs m))
(t (linalg-rank (grid:map-rows
(lambda (x)
(let ((c (grid:inner x zs)))
(map 'vector (lambda (a b) (- a (* c b))) x zs)))
vs)
m)))))))LINALG-RANKLet us test:
(time
(let ((n 120)
(m 150))
(print
(linalg-rank
(make-array
(list n m)
:initial-contents (loop repeat n collect
(loop repeat m collect
(- (random 1.0) 0.5))))))))120
Evaluation took:
1.252 seconds of real time
1.250837 seconds of total run time (1.214184 user, 0.036653 system)
[ Run times consist of 0.019 seconds GC time, and 1.232 seconds non-GC time. ]
99.92% CPU
87 lambdas converted
2,253,970,776 processor cycles
449,328,400 bytes consed