Add mpow function

[genmeblog]

M^n function implemented by mars0i

#data-science > Matrix-vector multiplication in fastmath

(defn mpow
  "Multiplies a square matrix by itself n times.  n must be a positive
  integer or positive floating-point number with no fractional part."
  [m n]
  (let [[h w] (fmat/shape m)]
    (cond
      (not= h w)
      (print (str "mpow only multiplies square matrices. Matrix shape was [" h " " w "].")) ; throw?

      (or (<= n 0)
          (not (zero? (rem n 1))))
      (print (str "mpow only accepts positive integer powers. Exponent was " n)) ; throw?

      :else
      (letfn [(either-pow [acc-mat k]
                (cond (= k 1) acc-mat
                      (even? k) (even-pow acc-mat k)
                      :else (fmat/mulm m (even-pow acc-mat (dec k)))))
              (even-pow [acc-mat k]
                (let [factor (either-pow acc-mat (/ k 2))]
                  (fmat/mulm factor factor)))]
        (either-pow m (long n)))))) ; allow float integers

[genmeblog]

As above works only for positive integers it's possible to make fractional power for non-singular matrices for which eigendecomposition exists $A = V D V^{-1}$, then $A^p = V D^p V^{-1}$, where $D^p$ contains a power of diagonal elements.