Multiple right-hand-sides

[jonasseglare]

When using core.matrix, it would be nice if we could solve linear systems and least-squares problems with multiple right hand sides. It currently does not seem to work, as demonstrated by my experiment. To reproduce it, do

lein new multiple-rhs

Then configure project file like this

(defproject multiple-rhs "0.1.0-SNAPSHOT"
  :description "FIXME: write description"
  :url "http://example.com/FIXME"
  :license {:name "Eclipse Public License"
            :url "http://www.eclipse.org/legal/epl-v10.html"}
  :dependencies [[org.clojure/clojure "1.6.0"]
                 [net.mikera/vectorz-clj "0.28.0"]
                 [net.mikera/core.matrix "0.32.1"]]
  :profiles {:dev {:plugins [[cider/cider-nrepl "0.8.1"]]}})

then call

lein deps

and write this in multiple_rhs/src/multiple_rhs/core.clj:

(ns multiple-rhs.core
  (:require [clojure.pprint :refer :all])
  (:refer-clojure :exclude [* - + == /])
  (:require [clojure.core.matrix :refer :all]
            [clojure.core.matrix.operators :refer :all])

  ; To get solve
  (:require [clojure.core.matrix.linear :refer :all]))

(set-current-implementation :vectorz)

; This is the left-hand-side of a linear system,
; a 2x2 matrix.
(def A [[-0.66718   3.05353] [0.60487   2.18721]])

; This is the right-hand-side of a linear system,
; a 2x3 matrix.
(def B [[-0.90260   0.49627  -0.74759]
        [1.02906  -0.75466   0.11171]])

; We would like (solve A B) to return this:
(mmul (inverse A) B)

; Ideally, we would like
(solve A B)
; which produces error:
; Cannot coerce INDArray with shape [I@4fa56922 to a vector

; We would like (least-squares A B) to return this:
; (mmul (inverse (mmul (transpose A) A)) (mmul (transpose A) B))

(least-squares A B)
; which produces error:
; Cannot coerce INDArray with shape [I@3364f8cf to a vector

and execute this code.

[jonasseglare]

By the way, here is the stack trace:

  Show: Clojure Java REPL Tooling Duplicates All  (24 frames hidden)

2. Unhandled clojure.lang.Compiler$CompilerException
   Error compiling:
   /home/jonas/programmering/clojure/multiple-rhs/src/multiple_rhs/core.clj:22:21

                 Compiler.java: 7142  clojure.lang.Compiler/load
                          REPL:    1  user/eval3331
                 Compiler.java: 6703  clojure.lang.Compiler/eval
                 Compiler.java: 6666  clojure.lang.Compiler/eval
                      core.clj: 2927  clojure.core/eval
                      main.clj:  239  clojure.main/repl/read-eval-print/fn
                      main.clj:  239  clojure.main/repl/read-eval-print
                      main.clj:  257  clojure.main/repl/fn
                      main.clj:  257  clojure.main/repl
                   RestFn.java: 1523  clojure.lang.RestFn/invoke
        interruptible_eval.clj:   67  clojure.tools.nrepl.middleware.interruptible-eval/evaluate/fn
                      AFn.java:  152  clojure.lang.AFn/applyToHelper
                      AFn.java:  144  clojure.lang.AFn/applyTo
                      core.clj:  624  clojure.core/apply
                      core.clj: 1862  clojure.core/with-bindings*
                   RestFn.java:  425  clojure.lang.RestFn/invoke
        interruptible_eval.clj:   51  clojure.tools.nrepl.middleware.interruptible-eval/evaluate
        interruptible_eval.clj:  183  clojure.tools.nrepl.middleware.interruptible-eval/interruptible-eval/fn/fn
        interruptible_eval.clj:  152  clojure.tools.nrepl.middleware.interruptible-eval/run-next/fn
                      AFn.java:   22  clojure.lang.AFn/run
       ThreadPoolExecutor.java: 1142  java.util.concurrent.ThreadPoolExecutor/runWorker
       ThreadPoolExecutor.java:  617  java.util.concurrent.ThreadPoolExecutor$Worker/run
                   Thread.java:  745  java.lang.Thread/run

1. Caused by java.lang.IllegalArgumentException
   Cannot coerce INDArray with shape [I@7610ea7c to a vector

                  Vectorz.java:  491  mikera.vectorz.Vectorz/toVector
                matrix_api.clj:  174  mikera.vectorz.matrix-api/avector-coerce*
                matrix_api.clj:  317  mikera.vectorz.matrix-api/eval13490/fn
                 protocols.clj:  847  clojure.core.matrix.protocols/eval8248/fn/G
                   default.clj: 1965  clojure.core.matrix.impl.default/eval10672/fn
                 protocols.clj:  847  clojure.core.matrix.protocols/eval8248/fn/G
                    linear.clj:  135  clojure.core.matrix.linear/solve
                      core.clj:   25  multiple-rhs.core/eval13210
                 Compiler.java: 6703  clojure.lang.Compiler/eval
                 Compiler.java: 7130  clojure.lang.Compiler/load
                          REPL:    1  user/eval3331
                 Compiler.java: 6703  clojure.lang.Compiler/eval
                 Compiler.java: 6666  clojure.lang.Compiler/eval
                      core.clj: 2927  clojure.core/eval
                      main.clj:  239  clojure.main/repl/read-eval-print/fn
                      main.clj:  239  clojure.main/repl/read-eval-print
                      main.clj:  257  clojure.main/repl/fn
                      main.clj:  257  clojure.main/repl
                   RestFn.java: 1523  clojure.lang.RestFn/invoke
        interruptible_eval.clj:   67  clojure.tools.nrepl.middleware.interruptible-eval/evaluate/fn
                      AFn.java:  152  clojure.lang.AFn/applyToHelper
                      AFn.java:  144  clojure.lang.AFn/applyTo
                      core.clj:  624  clojure.core/apply
                      core.clj: 1862  clojure.core/with-bindings*
                   RestFn.java:  425  clojure.lang.RestFn/invoke
        interruptible_eval.clj:   51  clojure.tools.nrepl.middleware.interruptible-eval/evaluate
        interruptible_eval.clj:  183  clojure.tools.nrepl.middleware.interruptible-eval/interruptible-eval/fn/fn
        interruptible_eval.clj:  152  clojure.tools.nrepl.middleware.interruptible-eval/run-next/fn
                      AFn.java:   22  clojure.lang.AFn/run
       ThreadPoolExecutor.java: 1142  java.util.concurrent.ThreadPoolExecutor/runWorker
       ThreadPoolExecutor.java:  617  java.util.concurrent.ThreadPoolExecutor$Worker/run
                   Thread.java:  745  java.lang.Thread/run

[mikera]

This looks like a very sensible enhancement to the current API.

Currently the core.matrix API doesn't require this support but I'm tempted to add it to vectorz-clj unilaterally :-)