mikera/vectorz-clj

Fast matrix and vector maths library for Clojure - as a core.matrix implementation

vectorz-clj

Fast vector and matrix library for Clojure, building on the Vectorz library and designed to work with the core.matrix array programming API.

vectorz-clj is designed so that you don't have to compromise, offering both:

• An idiomatic high-level Clojure API using core.matrix
• General purpose multi-dimensional arrays
• High performance (about as fast as you can get on the JVM). vectorz-clj is currently the fastest pure-JVM vector/matrix library available for Clojure

The library was originally designed for games, simulations and machine learning applications, but should be applicable for any situations where you need numerical double arrays.

Important features:

• "Pure" functions for an idiomatic functional programming style are provided. These return new vectors without mutating their arguments.
• Primitive-backed special purpose vectors and matrices for performance, e.g. Vector3 for fast 3D maths.
• Flexible DSL-style functions for manipulating vectors and matrices, e.g. the ability to create a "view" into a subspace of a large vector.
• core.matrix fully supported - see: https://github.com/mikera/core.matrix
• Pure cross-platform JVM code - no native dependencies
• "Impure" functions that mutate vectors are available for performance when you need it: i.e. you can use a nice functional style most of the time, but switch to mutation when you hit a bottleneck.

Documentation

vectorz-clj is intended to be used primarily as a core.matrix implementation. As such, the main API to understand is core.matrix itself. See the core.matrix wiki for more information:

• https://github.com/mikera/core.matrix/wiki

For more information about the specific details of vectorz-clj itself, see the vectorz-clj Wiki.

Status

vectorz-clj requires Clojure 1.4 or above, Java 1.7 or above, and an up to date version of core.matrix

vectorz-clj is reasonably stable, and implements all of the core.matrix API feature set.

License

Like Vectorz, vectorz-clj is licensed under the LGPL license:

• http://www.gnu.org/licenses/lgpl.html

Usage

Follow the instructions to install with Leiningen / Maven from Clojars:

• https://clojars.org/net.mikera/vectorz-clj

You can then use Vectorz as a standard core.matrix implementation. Example:

    (use 'clojure.core.matrix)
    (use 'clojure.core.matrix.operators)           ;; overrides *, + etc. for matrices
    
    (set-current-implementation :vectorz)  ;; use Vectorz as default matrix implementation
    
    ;; define a 2x2 Matrix
    (def M (matrix [[1 2] [3 4]]))
    M
    => #<Matrix22 [[1.0,2.0][3.0,4.0]]>
    
    ;; define a length 2 vector (a 1D matrix is considered equivalent to a vector in core.matrix)
    (def v (matrix [1 2]))
    v
    => #<Vector2 [1.0,2.0]>
    
    ;; Matrix x Vector elementwise multiply
    (mul M v)
    => #<Matrix22 [[1.0,4.0],[3.0,8.0]]>
    
    ;; Matrix x Vector matrix multiply (inner product)
    (inner-product M v)
    => #<Vector2 [5.0,11.0]>

For more examples see Wiki Examples

vectorz-clj's Issues

Which linear algebra algorithms need to implement

Hi, I saw your idea "Linear Algebra for Clojure" in wiki. And I want to make sure which algorithms need to implement? SVD, QR decomposition or others? Could you give a complete list of algorithms which need to implement? Thank you very much.

LUP decomposition seems to violate API

The API for LUP decomposition seems to indicate that

(defn LU-contract [A]
 (let [{:keys [L U P]} (lu A)]
  (= (mmul P A) (mmul L U))))

should return true.

When the implementation is vectorz, such as the case of

    [[ 1 -2  3]]
A =  [ 2 -5  12]
     [ 0  2 -10]]

it is not the case that (LU-contract A) returns true.

Instead, it appears to be the case that (= (mmul A) (mmul P L U)).

Not sure if this is a typo in the API docs or in the code itself - the tests in the tests folder seem to indicate that (= (mmul A) (mmul P L U)) is the intended behavior.

Failing case of QR decomposition

The following Qr decomposition case appears to fail:

[[  0 , 0 ],
 [ -1 , 0 ]]

with results

Q= [[  0.0,  1.0 ],
    [  1.0,  0.0 ]]
R= [[ -1.0,  0.0 ],
    [  0.0, -1.0 ]]

Worthwhile to add a project.clj?

I'd like to add a project.clj file so I can run the REPL and explore the code (as I look to modify it as discussed in the other issue). Would this be worthwhile, or would it conflict with other demands on the project?

If this isn't a good idea, then what is the recommended way to run a REPL against the code? (I don't know Maven very well.)

Operations on sparse matrices should produce sparse results

See core.matrix issue: mikera/core.matrix#115

set-indices returns nil

Using [net.mikera/vectorz-clj "0.45.0"]

(set-indices (matrix [[1 1] [2 2]]) [[1 1]] [9])
=> [[1 1] [2 9]]
(set-indices (matrix :vectorz [[1 1] [2 2]]) [[1 1]] [9])
=> nil

row-matrix for Vector returns nil instead of matrix

Found by Incanter's compliance test:

(bind-rows (matrix [1 2 3]) (matrix [4 5 6]))

Throws exception: java.lang.RuntimeException: Can't join an array to a nil value! because

incanter.core> (m/matrix [1 2 3])
#<Vector [1.0,2.0,3.0]>
incanter.core> (m/row-matrix (m/matrix [1 2 3]))
nil
incanter.core> (m/row-matrix [1 2 3])
#<Matrix [[1.0,2.0,3.0]]>

Works fine with other implementations: (clatrix, persistent-vector & NDArray:

incanter.core> (set-current-implementation :clatrix)
:clatrix
incanter.core> (bind-rows (matrix [1 2 3]) (matrix [4 5 6]))
 A 2x3 matrix
 -------------
 1.00e+00  2.00e+00  3.00e+00 
 4.00e+00  5.00e+00  6.00e+00 
incanter.core> (set-current-implementation :ndarray)
:ndarray
incanter.core> (bind-rows (matrix [1 2 3]) (matrix [4 5 6]))
#<NDArray [[1 2 3] [4 5 6]]>
incanter.core> (set-current-implementation :persistent-vector)
:persistent-vector
incanter.core> (bind-rows (matrix [1 2 3]) (matrix [4 5 6]))
[[1 2 3] [4 5 6]]

project.clj version ahead of pom.xml

In the develop branch, pom.xml says the version is 0.22.1-SNAPSHOT, while project.clj has 0.23.1-SNAPSHOT. I'm guessing that project.clj is wrong (did I do that?), since 0.22.0 is the latest release. Almost doesn't seem worth mentioning, but might as well.

(matrix/set-current-implementation :vectorz) or (set-current-implementation :vectorz); which one is correct?

What is the difference between both statements?

1. with a ns matrix, (matrix/set-current-implementation :vectorz)
2. without a ns matrix , (matrix/set-current-implementation :vectorz)

dependencies

[net.mikera/core.matrix "0.62.0"]
[net.mikera/vectorz-clj "0.48.0"]

Create a matrix successfully.

(use 'clojure.core.matrix)
(use 'clojure.core.matrix.random)
(use 'clojure.core.matrix.operators)  ;; for +, -, * etc.
**(matrix/set-current-implementation :vectorz)**
(def v (matrix [1 2]))

#'cljsl.ch1.examples/v

Create a matrix unsuccessfully.

(use 'clojure.core.matrix)
(use 'clojure.core.matrix.random)
(use 'clojure.core.matrix.operators)  ;; for +, -, * etc.
**(set-current-implementation :vectorz)**
(def v (matrix [1 2]))

Execution error at clojure.core.matrix.implementations/get-canonical-object (implementations.clj:64).
Unable to find implementation: [:vectorz]

Array creation on double array ~40x slower via core.matrix than when Vectorz/toVector used directly.

Here's a REPL transcript on the case of a 20-million long vector.
.
user> (def da (double-array (range 2e7)))

'user/da

;; 42 msec with direct call to Vectorz/toVector (Java API)
user> (time (def q (Vectorz/toVector da)))
"Elapsed time: 41.962 msecs"

'user/q

user> (time (def q (clojure.core.matrix/array da)))
Call to clojure.core.matrix.protocols/convert-to-nested-vectors took 1440.429 msec.
Call to mikera.vectorz.matrix-api/vectorz-coerce* took 1726.304 msec.
"Elapsed time: 1726.534 msecs"

'user/q

Most of the time seems to be spent converting the array into nested vectors (and, presumably, the rest is spent converting it back).

How easy or hard would it be to adapt this to avoid the round-trip in data types when a native double[](also, double[][] and higher dimensions) is passed in? It seems like this might require gut-level work on core.matrix in addition to vectorz-clj, but IMO it would be quite valuable to have fast vector/array creation when one already has a primitive double array. (I know those aren't Clojure idiomatic, but I find myself using them a lot for high-performance purposes.)

Bug in Cholesky

Cholesky decomposition seems to produce odd results where it should return nil, e.g.:

(cholesky (matrix :vectorz [[1 2] [2 1]]))
=> {:L* #<Matrix22 [[1.0,2.0],[0.0,NaN]]>, :L #<Matrix22 [[1.0,0.0],[2.0,NaN]]>}

"No implementation of method: :construct-matrix ..." on first load

When I run lein spec the matrix function fails with the error below. Running this as autotest lein spec -a it works on the secound run. So I think that something is wrong with my initialisation of the core.matrix implementation, because it seams to work on the second run, when the VM issn't restarted.

(ns athena.math-lib-assumptions-spec
  (:use [speclj.core]
        [clojure.core.matrix]
        [clojure.core.matrix.operators]))

(set-current-implementation :vectorz)

(describe "current-implementation"
  (it "is set as expected"
    (should= :vectorz (current-implementation))))


(describe "matrix"
  (with A (matrix [[1 2]
                   [3 4]]))

  (it "can be definded using standard clojure vector notation"
    (should= (matrix [[1 2]
                      [3 4]])
             @A))

1) matrix can be definded using standard clojure vector notation
     No implementation of method: :construct-matrix of protocol: #'clojure.core.matrix.protocols/PImplementation found for class: clojure.lang.Keyword
     java.lang.IllegalArgumentException: No implementation of method: :construct-matrix of protocol: #'clojure.core.matrix.protocols/PImplementation found for class: clojure.lang.Keyword

sinh with core.matrix

  (set-current-implementation :double-array)
  (pm (sinh (new-vector 10)))
  [0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000]

  (set-current-implementation :vectorz)
  ;; pm doesn't work with this
  (sinh (new-vector 10))
  #<Vector [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]>

set-nd on matrix returns nil

(mp/set-nd (matrix :vectorz [[1 2][3 4]]) [0 1] 3) => nil
Should fallback to calling its set-2d implementation
Come up in the mailing list https://groups.google.com/forum/#!topic/numerical-clojure/FsXo5j83hSM

Rationalize old API namespaces

The old vectorz-clj APIs do not match the preferred core.matrix API.

We should convert these to match core.matrix conventions.

Ideally these namespaces should contain exact equivalents of core.matrix functions, but avoid the protocol dispatch overhead (and use type hints where appropriate)

Create test.check generators

Should include test.check generators for different types of Vectorz arrays so that we can build a property-based testing suite.

This would be also useful as a general test of the core.matrix API

Missing support for PEigenDecomposition core.matrix protocol

It looks like Eigen and IEigenResult are included from the vectorz implementation, but the core.matrix.protocols.PEigenDecomposition protocol is not implemented. Any thoughts on an ETA? Thanks!

Infinity norm is incorrect

vectorz-clj/src/main/clojure/mikera/vectorz/matrix_api.clj

Lines 366 to 372 in 81be477

	(extend-protocol mp/PNorm
	INDArray
	(norm [m p]
	(cond
	(= java.lang.Double/POSITIVE_INFINITY p) (.elementMax m)
	(number? p) (Math/pow (.elementAbsPowSum m p) (/ 1.0 (double p)))
	:else (error "p must be a number"))))

This takes the .elementMax while it should be the max of the absolute values. (atm, the inf-norm can even be negative).
Unfortunately, .maxAbsElement this is not part of the INDArray interface. It's only implemented in AVector as far as I could tell.
A correct way that would work for INDArrays would be (.maxElement (.absCopy m)), but requires 2 traversals.

I guess the norm protocol can be implemented for INDArray and AVector, so it would be fast for vectors.

What is the best way to deal with this?

Element-wise division behavior is strange

Versions:

[net.mikera/core.matrix "0.13.1"]
[net.mikera/vectorz-clj "0.15.0"]

Faulty code:

linprog.core> ta
#<Matrix [[0.0,0.0,-1.0,-2.0],[1.0,-1.0,0.0,-1.0],[-1.0,0.0,-2.0,0.0]]>
linprog.core> (m/get-column ta 1)
#<StridedVector [0.0,-1.0,0.0]>
linprog.core> tb
#<Vector3 [1.0,3.0,0.0]>
linprog.core> (m/div tb (m/get-column ta 1))
#<Vector3 [1.0,3.0,0.0]>
linprog.core> (m/div tb [0.0,-1.0,0.0])
#<Vector3 [1.0,3.0,0.0]>
linprog.core> (m/div [1.0,3.0,0.0] [0.0,-1.0,0.0])
ArithmeticException Divide by zero  clojure.lang.Numbers.divide (Numbers.java:156)

matrix->rows->matrix round trip does not work

After the update to 0.45.0, the following code throws an exception:

(def m (matrix/matrix :vectorz [[1 2 3][4 5 6]])
(def rows (matrix/rows m))
(def m2 (matrix/matrix :vectorz (vec rows)))

The exception thrown is: java.lang.UnsupportedOperationException count not supported on this type: ArraySubVector. A similar exception with StridedVector is thrown when I use matrix/columns instead of matrix/rows.

I don't see this problem with 0.44.1 and if I remove :vectorz from one or both of the above matrix/matrix calls then it works.

lt and friends return bogus results

I'm getting the following:

> (matrix/lt (matrix/array [1 2 -1 2]) 0.0)
#vectorz/vector [0.0,0.0,0.0,0.0] ;; I'd expect the third element to be 1

> (matrix/lt [1 2 -1 2] 0)
[0 0 1 0] ;; this is okay

> (matrix/gt (matrix/array [1 2 -1 2]) 0.0)
#vectorz/vector [1.0,1.0,1.0,1.0] ;; at least it's consistent with gt...

What gives? I'm on vectorz-clj 0.44.1.

pm (pretty-print function): error on vectorz matrices

=> (pm (matrix :persistent-vector [[1 2][3 4]]))
[[1.000 2.000]
 [3.000 4.000]]
nil
=> (pm (matrix :ndarray [[1 2][3 4]]))
[[1.000 2.000]
 [3.000 4.000]]
nil
=> (pm (matrix :clatrix [[1 2][3 4]]))
[[1.000 2.000]
 [3.000 4.000]]
nil
=> (pm (matrix :vectorz [[1 2][3 4]]))
ClassCastException mikera.vectorz.impl.VectorIndexScalar cannot be cast to java.lang.Number  clojure.lang.RT.doubleCast (RT.java:1222)

basic bug with sparse matrix operations

cn.product.index> (matrix/sparse-array [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] )
ArrayIndexOutOfBoundsException 419 mikera.vectorz.Vector.getElements (Vector.java:194)
cn.product.index> (matrix/sparse-matrix [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] )
ArrayIndexOutOfBoundsException 419 mikera.vectorz.Vector.getElements (Vector.java:194)
cn.product.index>

Support higher dimension selects

Currently the PSelect protocol is only implemented for 1D vectors.

We should support the PSelect protocol on arrays of any size.

Overridden operators?

I get this warning message:

"WARNING: == already refers to: #'clojure.core/== in namespace: music.core, being replaced by: #'clojure.core.matrix.operators/==\nWARNING: + already refers to: #'clojure.core/+ in namespace: music.core, being replaced by: #'clojure.core.matrix.operators/+\nWARNING: * already refers to: #'clojure.core/* in namespace: music.core, being replaced by: #'clojure.core.matrix.operators/*\nWARNING: - already refers to: #'clojure.core/- in namespace: music.core, being replaced by: #'clojure.core.matrix.operators/-\nWARNING: / already refers to: #'clojure.core// in namespace: music.core, being replaced by: #'clojure.core.matrix.operators//\n

Vectorz breaks dot product for matricies

Hi there,

In vectorz-clj 0.46.0 you can't use the matrix/dot function on vectorz' matrices but this works in vanilla core.matrix.

Core.Matrix

> (use 'clojure.core.matrix)
;;=> nil
> (def v (array [[1 2 3] [4 5 6]]))
;;=> #`user/v
> (dot v v)
;;=> [[9 12 15] [24 33 42]]

Vectorz

> (set-current-implementation :vectorz)
;;=> :vectorz
> (def v (array [[1 2 3] [4 5 6]]))
;;=> #'user/v
> (dot v v)
;;=> RuntimeException Can't coerce to AVector with shape: [2 3]  mikera.vectorz.matrix-api/avector-coerce* (matrix_api.clj:198)

Add license to readme

From reading the pom, it seems that it's licensed under LGPL.
In any case, it would be really nice if the license was mentioned
in the readme, as it's not really possible to use a library without
knowing what the license is.

e= fails to detect unequality

I stumbled upon this behavior:

e= seams failing to detect if two matrices are NOT equal, while the std = function does.

Is this indeed a bug or am I wrong about my assumptions?

(ns athena.math-lib-assumptions-spec
  (:use [speclj.core]
        [clojure.core.matrix]
        [clojure.core.matrix.operators]))

(set-current-implementation :vectorz)


(describe "e="
  (with expected (matrix [[1 2]
                          [3 4]]))

  (it "detects equality"
    (should (e= @expected
                (matrix [[1 2]
                         [3 4]]))))
  (it "FAILS TO DETECTS UNEQUALITY DONT USE IT"
    (should (e= (matrix @expected)
                (matrix [[2 2]
                         [3 4]])))))


(describe "= for matrices"
  (with expected (matrix [[1 2]
                          [3 4]]))

  (it "detects equality"
    (should= @expected
              (matrix [[1 2]
                       [3 4]])))
  (it "detects unequality"
    (should-not= (matrix @expected)
                 (matrix [[2 2]
                          [3 4]]))))

emin returns 1.7976931348623157E308 when minimum is Infinity

The emin function in 0.48.0 returns 1.7976931348623157E308 (max double) when the minimum of a matrix is, in fact, Infinity.

Quick example:

(-> [##Inf] matrix/matrix matrix/emin)

returns 1.7976931348623157E308.

cast mikera.vectorz.Vector to clojure.lang.PersistentVector

Is there anyway to do the above typecast ?

I am trying to use incanter.charts/scatter-plot

The following works fine:

(incanter.charts/scatter-plot [1 2 3] [4 5 6])

But this:

(require [clojure.core.matrix :as m])
(m/set-current-implementation :vectorz)

(def x (m/matrix (for [x (range 10)] (+ x (* 0.1 (rand))))))
(def y (m/matrix (for [x (range 10)] (+ x (* 0.1 (rand))))))
(incanter.charts/scatter-plot x y)

throws the following error:

Exception in thread "main" java.lang.ClassCastException: mikera.vectorz.Vector cannot be cast to java.lang.Number, compiling:(/tmp/form-init6889720767473931280.clj:1:72)
        at clojure.lang.Compiler.load(Compiler.java:7142)
        at clojure.lang.Compiler.loadFile(Compiler.java:7086)
        at clojure.main$load_script.invoke(main.clj:274)
        at clojure.main$init_opt.invoke(main.clj:279)
        at clojure.main$initialize.invoke(main.clj:307)
        at clojure.main$null_opt.invoke(main.clj:342)
        at clojure.main$main.doInvoke(main.clj:420)
        at clojure.lang.RestFn.invoke(RestFn.java:421)
        at clojure.lang.Var.invoke(Var.java:383)
        at clojure.lang.AFn.applyToHelper(AFn.java:156)
        at clojure.lang.Var.applyTo(Var.java:700)
        at clojure.main.main(main.java:37)
Caused by: java.lang.ClassCastException: mikera.vectorz.Vector cannot be cast to java.lang.Number
        at clojure.lang.RT.doubleCast(RT.java:1222)
        at incanter.charts$scatter_plot_STAR_$fn__13382.invoke(charts.clj:1365)
        at clojure.core$map$fn__4249.invoke(core.clj:2564)
        at clojure.lang.LazySeq.sval(LazySeq.java:40)
        at clojure.lang.LazySeq.seq(LazySeq.java:49)
        at clojure.lang.RT.seq(RT.java:484)
        at clojure.core$seq.invoke(core.clj:133)
        at clojure.core$dorun.invoke(core.clj:2855)
        at incanter.charts$scatter_plot_STAR_.doInvoke(charts.clj:1363)
        at clojure.lang.RestFn.applyTo(RestFn.java:142)
        at clojure.core$apply.invoke(core.clj:624)
        at clj_lin_reg.core$_main.invoke(core.clj:13)
        at clojure.lang.Var.invoke(Var.java:375)
        at user$eval5.invoke(form-init6889720767473931280.clj:1)
        at clojure.lang.Compiler.eval(Compiler.java:6703)
        at clojure.lang.Compiler.eval(Compiler.java:6693)
        at clojure.lang.Compiler.load(Compiler.java:7130)
        ... 11 more

What I had in mind was to cast it to a persistent vector before transforming it. Is there an API for this library to do that?

orz

symmetric? with INDArray ?

In order to add a symmetric? predicate to PMatrixPredicates in matrix_api.clj, Mike has pointed out that I can use .isSymmetric in AMatrix, which works well.

The extend-protocol for PMatrixPredicates also has a condition for INDArray. INDArray doesn't have isSymmetric(), however. Should it? Under what conditions is a vectorz-clj matrix INDArray but not AMatrix? I don't know how to test this for this case. Should I just use a pure Clojure test in this condition?

(Is this the place to ask these questions? Not sure if this counts as an issue per se.)

Issues and inconsistencies in making block matrices

I was forced to switch implementations b/c for some reason the default implementation doesn't implement the full matrix.core API (specifically the lu function)

However once I switched to vectorz my old code stopped working

> (set-current-implementation :persistent-vector)
:persistent-vector
> (block-diagonal-matrix [(identity-matrix 2) (identity-matrix 2)])
[[1.0 0.0 0.0 0.0]
 [0.0 1.0 0.0 0.0]
 [0.0 0.0 1.0 0.0]
 [0.0 0.0 0.0 1.0]]
> (set-current-implementation :vectorz)
:vectorz
> (block-diagonal-matrix [(identity-matrix 2) (identity-matrix 2)])
Execution error (ClassCastException) at clojure.core.matrix.impl.defaults/eval16737$fn$aux$iter$fn$fn (defaults.cljc:2202).
class mikera.matrixx.impl.IdentityMatrix cannot be cast to class clojure.lang.IFn (mikera.matrixx.impl.IdentityMatrix and clojure.lang.IFn are in unnamed module of loader 'app')

non-zero-indices fails on >1d arrays

The default implementation (in default.clj) of PNonZeroIndices will invoke the array/matrix as a function when attempting to get the non-zero indices for an inner dimension.

Successful broadcasting by comparison operators depends on operand order

(require '[clojure.core.matrix :as m])

(m/set-current-implementation :vectorz)

(def vect (m/array [1 0]))
(def mat (m/array [[0 0] [1 1]]))

(m/eq mat vect) ;; => [[0 1] [1 0]]
(m/eq vect mat) ;; java.lang.IllegalArgumentException
                ;; Can't broadcast Matrix with shape [2,2] to shape: [2]

The above exception does not occur with the ndarray core.matrix implementation. With vectorz, it occurs with at least some of the other comparison operators (e.g. ne, gt) but not with at least some of the other arithmetic operators (e.g. add).

QR fails with rows < columns

(qr (matrix :vectorz [[1 2]]))
IllegalArgumentException Wrong matrix size: rows < columns  mikera.matrixx.decompose.QR.decomposeInternal (QR.java:66)

I think this should succeed with something like the following result:
Q= [[1.0]]
R= [[1.0 2.0]]

This seems to be a perfectly valid QR decomposition?

Multiple right-hand-sides

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.

Documentation error: clojure.core.matrix.operators/* implements element-wise multiplication.

Current documentation indicates that (* M v), where M is a matrix and v is a vector, will perform linear-algebraic matrix multiplication. In fact, it performs an element-wise multiplication, like so:

user=> (use 'clojure.core.matrix)
nil
user=> (use 'clojure.core.matrix.operators)
WARNING: / already refers to: #'clojure.core// in namespace: user, being replaced by: #'clojure.core.matrix.operators//
WARNING: - already refers to: #'clojure.core/- in namespace: user, being replaced by: #'clojure.core.matrix.operators/-
WARNING: * already refers to: #'clojure.core/* in namespace: user, being replaced by: #'clojure.core.matrix.operators/*
WARNING: + already refers to: #'clojure.core/+ in namespace: user, being replaced by: #'clojure.core.matrix.operators/+
WARNING: == already refers to: #'clojure.core/== in namespace: user, being replaced by: #'clojure.core.matrix.operators/==
nil
user=> (set-current-implementation :vectorz)
:vectorz
user=> (def M (matrix [[1 2] [3 4]]))
#'user/M
user=> (def v (matrix [1 2]))
#'user/v
user=> (* M v)
#<Matrix22 [[1.0,4.0],[3.0,8.0]]>

join-along throws exception

The following code works for all built-in implementations of core.matrix:

(require '[clojure.core.matrix :as m])
(let [a (m/matrix [[1 3 5] [2 4 6]])
      b (m/matrix [7 8])]
  (m/join-along 1 a b))
=> [[1 3 5 7]
    [2 4 6 8]]

However, with the :vectorz implementation, I get the following exception:

IllegalArgumentException Incompatible shapes:
Matrix with shape [2,3] vs. Vector with shape [2] 
mikera.arrayz.impl.JoinedArray.join (JoinedArray.java:33)

IllegalStateExceptions may be caused by AOT'ed code in Clojure 1.11

Found this bug and (somewhat) a solution, posting here for reference and in case this causes an issue for someone down the road.

The error, from calling clojure.core.matrix.linear/norm is:
java.lang.IllegalStateException: Attempting to call unbound fn: #'clojure.core.matrix.impl.mathsops/abs

Full Stacktrace and System Info

java.lang.IllegalStateException: Attempting to call unbound fn: #'clojure.core.matrix.impl.mathsops/abs
        at clojure.lang.Var$Unbound.throwArity(Var.java:45)
        at clojure.lang.AFn.invoke(AFn.java:32)
        at clojure.core$map$fn__5935.invoke(core.clj:2770)
        at clojure.lang.LazySeq.sval(LazySeq.java:42)
        at clojure.lang.LazySeq.seq(LazySeq.java:51)
        at clojure.lang.LazySeq.first(LazySeq.java:73)
        at clojure.lang.RT.first(RT.java:692)
        at clojure.core$first__5449.invokeStatic(core.clj:55)
        at clojure.core$first__5449.invoke(core.clj:55)
        at clojure.core.matrix.impl.sequence$fn__8123.invokeStatic(sequence.cljc:117)
        at clojure.core.matrix.impl.sequence$fn__8123.invoke(sequence.cljc:114)
        at clojure.core.matrix.protocols$fn__468$G__420__473.invoke(protocols.cljc:75)
        at clojure.core.matrix.impl.persistent_vector$persistent_vector_coerce.invokeStatic(persistent_vector.cljc:96)
        at clojure.core.matrix.impl.persistent_vector$persistent_vector_coerce.invoke(persistent_vector.cljc:93)
        at clojure.core.matrix.impl.persistent_vector$fn__5183.invokeStatic(persistent_vector.cljc:138)
        at clojure.core.matrix.impl.persistent_vector$fn__5183.invoke(persistent_vector.cljc:125)
        at clojure.core.matrix.protocols$fn__381$G__318__388.invoke(protocols.cljc:40)
        at clojure.core.matrix.impl.common$construct_matrix.invokeStatic(common.cljc:52)
        at clojure.core.matrix.impl.common$construct_matrix.invoke(common.cljc:48)
        at clojure.core.matrix.impl.common$mapmatrix.invokeStatic(common.cljc:70)
        at clojure.core.matrix.impl.common$mapmatrix.invoke(common.cljc:58)
        at clojure.core.matrix.impl.persistent_vector$fn__5735.invokeStatic(persistent_vector.cljc:651)
        at clojure.core.matrix.impl.persistent_vector$fn__5735.invoke(persistent_vector.cljc:637)
        at clojure.core.matrix.protocols$fn__3932$G__3923__3953.invoke(protocols.cljc:1041)
        at clojure.core.matrix.impl.defaults$fn__8013.invokeStatic(defaults.cljc:2418)
        at clojure.core.matrix.impl.defaults$fn__8013.invoke(defaults.cljc:2413)
        at clojure.core.matrix.protocols$fn__4704$G__4699__4711.invoke(protocols.cljc:1225)
        at clojure.core.matrix.linear$norm.invokeStatic(linear.cljc:21)
        at clojure.core.matrix.linear$norm.invoke(linear.cljc:9)
        at clojure.core.matrix.linear$norm.invokeStatic(linear.cljc:20)
        at clojure.core.matrix.linear$norm.invoke(linear.cljc:9)
		...

Tested on previously functioning code with

• Linux
• OpenJDK Runtime Environment Temurin-17.0.5+8
• Clojure 1.11.1

---

This is a known bug (CLJ-2712) and fixed in Clojure 1.12.
I was able to use the 1.12.0-alpha1 release and it works fine now.

I have not tried any AOT solutions, nor other versions of Clojure, nor other JVMs.

Inconsistency with QR results

The default implementation and vectorz-clj implementations of QR decomposition produce quite different shapes of result for rectangular matrices:

With Clojure persistent vectors:

(qr [[1 2] [3 4] [5 6] [7 8]])
  {:R [[-9.165151389911681 -10.910894511799619] 
       [0 -0.975900072948533] 
       [0 0] 
       [0 0]], 
    :Q [[-0.10910894511799607 -0.8295150620062532 -0.3945010222838296 -0.3799591338775966] 
        [-0.32732683535398854 -0.4391550328268402 0.24279654570435638 0.800655879510063] 
        [-0.5455447255899809 -0.04879500364742667 0.6979099754427761 -0.46143435738733596] 
        [-0.7637626158259733 0.34156502553198664 -0.5462054988633029 0.040737611754869674]]}

With vectorz-clj:

(qr (matrix :vectorz [[1 2] [3 4] [5 6] [7 8]]))
{:R #<Matrix [[-9.16515138991168,-10.91089451179962],
              [0.0,-0.9759000729485332]]>, 
 :Q #<Matrix [[-0.10910894511799629,-0.8295150620062537],
              [-0.3273268353539886,-0.43915503282683965],
              [-0.545544725589981,-0.048795003647426116],
              [-0.7637626158259734,0.3415650255319861]]>}

I think the API should specify what shape of result should be returned, or if there is a choice of result format these should be specified through the decomposition options.

Exception creating a sparse array

Hi,
I'm running into an ArrayIndexOutOfBoundsException when trying to create a sparse-array.
These are the versions I use (latest stable at this moment):
[net.mikera/core.matrix "0.52.0"]
[net.mikera/vectorz-clj "0.44.0"]

When I read the edn data from the attached zip file, and try to create a sparse-array, I get an exception:

> (def fd (read-string (slurp "failing-data.edn")))
> (m/sparse-array fd)
ArrayIndexOutOfBoundsException 8498  mikera.vectorz.impl.AStridedVector.getElements (AStridedVector.java:306)

It's not due to the matrix not being rectangular (although even non-rectangular seems to be supported):

> (distinct (map count fd))
(128)

Are there some restrictions on the creation of sparse-arrays? I didn't see any in the docs.
If not, I assume this is a vectorz-clj issue, rather than a core.matrix issue.

Thanks,
-Mathias

failing-data.zip

The inverse function in combination with the diagonal-matrix function seems to be bugged

If I run
(inverse (matrix [[1 0] [0 2]]))
I get
#vectorz/matrix [[1.0,-0.0],[-0.0,0.5]]

However when I run
(inverse (diagonal-matrix [1 2]))
I get
#vectorz/matrix [[Infinity,0.0],[0.0,Infinity]]

Not sure what is going on here because the output of both
(matrix [[1 0] [0 2]]) and (diagonal-matrix [1 2]) is #vectorz/matrix [[1.0,0.0],[0.0,2.0]]

function nth not supported

Hi,

Awesome library. I think I found a bug however.

Neither type Matrix nor Vector seem to support the function nth. This has the result that the type needs to be converted – e.g. using to-nested-vectors – before destructuring can be used.

Examples:

(nth (m/array [1 1]) 0)

and

(let [[x y] (m/matrix [2 2])]
    x)

both return the error message

Exception while evaluating the expression.
#unrepl/browsable
  [#error {:via [{:type java.lang.UnsupportedOperationException,
                  :message "nth not supported on this type: Vector",
                  :at [clojure.lang.RT nthFrom "RT.java" 928]}],
           :trace [[clojure.lang.RT nthFrom "RT.java" 928]
                   [clojure.lang.RT nth "RT.java" 893]
                   [user$eval34322 invokeStatic "NO_SOURCE_PATH" 1]
                   [user$eval34322 invoke "NO_SOURCE_PATH" 1]
                   [clojure.lang.Compiler eval "Compiler.java" 7176]
                   [clojure.lang.Compiler eval "Compiler.java" 7131]
                   [clojure.core$eval invokeStatic "core.clj" 3214]
                   [clojure.core$eval invoke "core.clj" 3210]
                   [user$eval34315 invokeStatic "NO_SOURCE_PATH" 52]
                   [user$eval34315 invoke "NO_SOURCE_PATH" 52] /117],
           :cause "nth not supported on this type: Vector"}
   /118]

Is this intended?

Cheers

V* in SVD should be transposed

I am not 100% sure if vectorz does it wrong or other libs, I am continuing to investigate this. However, here is the issue, illustrated with an example:

user=> (clojure.core.matrix/set-current-implementation :vectorz)
:vectorz
(:V* (decomp-svd (covariance [[1 2 3 ] [ 4 5 6 ] [7 8 9] [10 11 12]])))
#<mikera.matrixx.Matrix@492d4c90
[[0.5773502691896257,0.5773502691896255,0.5773502691896257],
[0.816496580927726,-0.40824829046386285,-0.40824829046386296],
[0.0,-0.7071067811865475,0.7071067811865475]]>

user=> (clojure.core.matrix/set-current-implementation :clatrix)
:clatrix
user=> (:V* (decomp-svd (covariance [[1 2 3 ] [ 4 5 6 ] [7 8 9] [10 11 12]])))
((-0.5773502691896256 0.8164965809277263 0.0)
 (-0.5773502691896258 -0.40824829046386285 -0.7071067811865476)
 (-0.577350269189626 -0.40824829046386296 0.7071067811865475))

This is done via Incanter / core.matrix. As you can see the resulting V* are different as one is the transpose of the other. There are also differences in the signs, but they are not significant for SVD and we don't need to bother.

R returns the same data as the clatrix implementation:

> svd(cov(dat))
...
$v
           [,1]       [,2]       [,3]
[1,] -0.5773503  0.0000000  0.8164966
[2,] -0.5773503 -0.7071068 -0.4082483
[3,] -0.5773503  0.7071068 -0.4082483

Thus I think vectorz does it wrong. The other parts of the return value (U,S) are correct (aka the same) in all implementations and R.

Exponential with matrices

Would it be possible to add support for doing exponential with matrices.

In python numpy, it is possible to do

$ A = np.array([[1,2],[3,4]])
$ 2 ** A
array([[ 2,  4],
       [ 8, 16]])

# as well as
$ A ** 2
array([[ 1,  4],
       [ 9, 16]])

This would be super-useful in more complex equations. I'm looking forward to hear from you!

Issue initialising large matrices with data

Using:

[net.mikera/core.matrix "0.29.1"]
[net.mikera/vectorz-clj "0.25.0"]
[org.clojure/clojure "1.6.0"]

Attempting to initialise a "large" matrix does not preserve the provided values (all cells are set to 0.0)

(require '[clojure.core.matrix :as mat]) ;; => nil
(mat/current-implementation) ;; => :persistent-vector
(let [w 200 h 500] (mat/esum (mat/matrix (for [y (range h)] (for [x (range w)] (+ x (* y w))))))) ;; => 4999950000
(let [w 250 h 500] (mat/esum (mat/matrix (for [y (range h)] (for [x (range w)] (+ x (* y w))))))) ;; => 7812437500

(mat/set-current-implementation :vectorz) ;; => :vectorz
(let [w 200 h 500] (mat/esum (mat/matrix (for [y (range h)] (for [x (range w)] (+ x (* y w))))))) ;; => 4.99995E9
(let [w 250 h 500] (mat/esum (mat/matrix (for [y (range h)] (for [x (range w)] (+ x (* y w))))))) ;; => 0.0

set-column throws IndexOutOfBoundsException if index = column-count

In vectorz-clj 0.45.0, trying to set a column in column n for an array of size m*n throws an IndexOutOfBoundsException.
Sample code:

user> (matrix/set-current-implementation :vectorz)
:vectorz
user> (-> (matrix/new-array [2 3])
          (matrix/set-column 2 [3 4]))
IndexOutOfBoundsException Matrix with shape [2,3] does not have slice: 2 on dimension 0  mikera.matrixx.impl.ARectangularMatrix.checkRow (ARectangularMatrix.java:93)

The exact same code works as expected with clatrix:

user> (-> (matrix/new-array [2 3])
          (matrix/set-column 2 [3 4]))
 A 2x3 matrix
 -------------
 0,00e+00  0,00e+00  3,00e+00 
 0,00e+00  0,00e+00  4,00e+00 

Setting a column also works without problems for positions 0<n.