The Kitchen Sink and Other Oddities

Atabey Kaygun

Calculating The Correct Rank of a Matrix

Description of the problem

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.

Linear algebraic rank of a matrix

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-RANK

Let 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