fast linear algebra for Java http://jblas.org Mikio L. Braun
[email protected] TU Berlin, Franklinstr. 28/29, 10587 Berlin
June 23, 2010
Overview Main features:
Matrix library for Java based on native BLAS and LAPACK.
Support for single and double precision floats, real and complex matrices.
Vectors and two-dimensional dense matrices only.
For large matrix practically native performance (for example, 10 GFLOPS on a 2GHz Intel Core2Duo processor for matrix-matrix multiplication).
Precompiled “fat jar” for Linux (Intel), Mac OS, and Windows.
Parses FORTRAN code to automatically generate JNI stubs.
Benchmarks
Structure Main classes
FloatMatrix, DoubleMatrix – real matrices.
ComplexFloatMatrix, ComplexDoubleMatrix – complex matrices.
Computational routines as static methods
Eigen – eigendecomposition
Solve – solving linear equations
Singular – singular value decomposition
Decompose – decompositions like LU, Cholesky, ...
Geometry – centering, normalizing, ...
Background on Native Calls in Java
Native calls are expensive:
Arrays are always copied on native calls.
For large arrays lots of cache misses.
Doesn’t make sense to call native code if operation is O (1)! (1)!
What about direct buffers?
Direct buffers aren’t garbage collected properly!
Only makes sense to have a fixed number of direct buffers used for copying → directly use JNI with arrays.
Nevertheless, Java is as fast as C for simple things like vector addition!
Matrix creation
Doubl DoubleMa eMatr trix ix a = new Doubl DoubleMa eMatr trix ix(10 (10, , 5); // 10 * 5 matrix Doubl DoubleMa eMatr trix ix x = new Doubl DoubleMa eMatr trix ix(10 (10); ); // vec vector of leng ength 10 Doubl DoubleMa eMatr trix ix y = Doubl DoubleMa eMatr trix. ix.ze zeros ros(1 (10, 0, 5); 5); Doubl DoubleMa eMatr trix ix z = Doubl DoubleMa eMatr trix. ix.on ones( es(3, 3, 4); Doubl DoubleMa eMatr trix ix g = Doubl DoubleMa eMatr trix. ix.ra randn ndn(3 (3, , 4); 4); DoubleMa DoubleMatri trix x u = DoubleMa DoubleMatrix trix.ra .rand(3 nd(3); );
Dimensions: rows, columns, length
Duplicating and copying: dup(), copy()
Transposing: transpose()
Element access Doubl DoubleMa eMatr trix ix a = new Doubl DoubleMa eMatr trix ix(10 (10, , 5); // acce access ssin ing g elem elemen ents ts by row row and and colu column mn a.put a.put(3, (3, 2, 10.0) 10.0); ; a.get a.get(2, (2, 3); 3); // acce access ssin ing g elem elemen ent t in line linear ar orde order r a.get(15); a.put(20 a.put(20, , 1.0); 1.0); // for for exam example ple, , for imple impleme menti nting ng eleme elementw ntwis ise e opera operati tion ons: s: for (int i = 0; i < a.length; i++) a.pu a.put( t(i, i, a.ge a.get( t(i) i) * 3); 3);
Accessing rows and columns: putRow(), getRow(), ...
Pass a buffer object to prevent spurious object creation: Doubl DoubleMa eMatr trix ix a = Doubl DoubleMa eMatr trix. ix.ra randn ndn(1 (10, 0, 100); 100); Doubl DoubleMa eMatr trix ix buf buf = new new Doub DoubleM leMat atrix rix(1 (10); 0); for (int i = 0; i < a.columns; i++) { a.getColumn(i, a.getColumn(i, buf); // do some someth thin ing g with with buf buf }
Simple Arithmetic Doub Double leMa Matr trix ix a = new new Doub Double leMa Matr trix ix(3 (3, , 3, 1, 2, 3, 4, 5, 6, 7, 8, 9); Doub Double leMa Matr trix ix x = new new Doub Double leMa Matr trix ix(3 (3, , 1, 10, 10, 11, 11, 12); 12); Doubl DoubleMa eMatr trix ix y = a.mmu a.mmul(x l(x); ); Doubl DoubleMa eMatr trix ix z = x.add x.add(y) (y); ;
Supported Operations:
Basic arithmetics: add, sub, mul, mmul, div, dot
Element-wise logical operations: and, or, xor, not, lt, le, gt, ge, eq, ne
Rows and vectors to a matrix: addRowVector, addColumnVector, ...
Machine Learning: Kernel Ridge Regression with Gaussian Kernel Let’s implement the following: Noisy sinc data set learned with Kernel Ridge Regression with a Gaussian Kernel. Sinc Function unction
sinc(x ) =
sin(x )/x 1
x = 0 x = 0.
DoubleMa DoubleMatri trix x sinc(Do sinc(Double ubleMat Matrix rix x) { return sin(x).div(x); sin(x).div(x); }
Not safe, what about x = 0? DoubleMa DoubleMatri trix x safeSin safeSinc(Do c(Doubl ubleMat eMatrix rix x) { return sin(x).div(x.ad sin(x).div(x.add(x.eq( d(x.eq(0))); 0))); }
Computing a sinc data set
X ∼ uniformly from -4, 4 Y = sinc(x ) + σ 2 ε ε
ε ∼ N (0, (0, 1). 1). Doubl DoubleMa eMatr trix ix[] [] sincD sincDat atase aset( t(int int n, doub double le noise noise) ) { DoubleM DoubleMatr atrix ix X = rand(n) rand(n).mul .mul(8) (8).sub .sub(4); (4); Doub DoubleM leMat atrix rix Y = sinc( sinc(X) X) .add( .add( randn randn(n) (n).m .mul( ul(no noise ise) ) ); retu return rn new new Doub DoubleM leMat atrix rix[] [] {X, {X, Y}; }
Gaussian kernel
k (x , z ) = exp −
2
x − z w
DoubleMa DoubleMatri trix x gaussia gaussianKer nKernel nel(dou (double ble w, DoubleM DoubleMatr atrix ix X, Doub DoubleM leMat atrix rix Z) { Doub Double leMa Matr trix ix d = Geometry.pairwiseSquaredDistances(X.transpose(), Z.transpose()); return exp(d.div(w).ne exp(d.div(w).neg()); g()); }
Kernel Ridge Regression
KRR learns a “normal” kernel model n
f (x ) =
k (x , X i i )αi
i =1
with α = (K + λI )
−1
Y ,
where K is the kernel matrix K ij ij = k (X i i , X j ).
DoubleMa DoubleMatri trix x learnKR learnKRR(Do R(Doubl ubleMat eMatrix rix X, DoubleMa DoubleMatrix trix Y, doub double le w, doub double le lamb lambda da) ) { int n = X.row rows; Doub DoubleM leMat atrix rix K = gauss gaussia ianKe nKern rnel el(w, (w, X, X); X); K.addi(eye(n).muli(lambda)); DoubleM DoubleMatr atrix ix alpha alpha = Solve.s Solve.solv olveSym eSymmetr metric(K ic(K, , Y); return return alpha; alpha; } DoubleMa DoubleMatri trix x predict predictKRR( KRR(Dou DoubleM bleMatri atrix x XE, DoubleM DoubleMatri atrix x X, doub double le w, Doub Double leMa Matr trix ix alph alpha) a) { Doub DoubleM leMat atrix rix K = gauss gaussia ianKe nKern rnel el(w, (w, XE, X); X); return K.mmul(alpha); K.mmul(alpha); }
Computing the mean-squared error
1 n
n
(Y i i − Y i )2
i =1
doubl double e mse( mse(Dou Doubl bleMa eMatr trix ix Y1, Y1, Doub DoubleM leMat atrix rix Y2) Y2) { Doub DoubleM leMat atrix rix diff diff = Y1.su Y1.sub(Y b(Y2) 2); ; return return pow(dif pow(diff, f, 2).mean 2).mean(); (); }
Conjugate Gradients 1. r ← b − Ax 2. p ← r 3. repeat r T r p T Ap
4.
α←
5.
x ← x + αp
6.
r ← r − αAp
7.
if r is sufficiently small, exit loop
T
r r r T r
8.
β ←
9.
p ← r + β p p
10.
r ← r .
11. end repeat
Double DoubleMat Matrix rix cg(Dou cg(Double bleMat Matrix rix A, Double DoubleMat Matrix rix b, Doub Double leMa Matr trix ix x, doub double le thre thresh sh) ) { int int n = x.le x.leng ngth th; ; Doub Double leMa Matr trix ix r = b.su b.sub( b(A. A.mm mmul ul(x (x)) )); ; // 1 DoubleMatrix p = r.dup(); // 2 double alpha = 0, beta = 0; Doub Double leMa Matr trix ix r2 = zero zeros( s(n) n), , Ap = zero zeros( s(n) n); ; while (true) { // 3 A.mmul A.mmuli(p i(p, , Ap); Ap); alpha = r.dot(r) / p.dot(Ap); // 4 x.addi(p.mul(alpha)); // 5 r.subi(Ap.mul(alpha), r2); // 6 double error = r2.norm2(); // 7 System System.ou .out.p t.prin rintf( tf("Re "Resid sidual ual error error = %f\n", %f\n", error); error); if (err (error or < thre thresh sh) ) break; beta = r2.dot(r2) / r.dot(r); // 8 r2.addi(p.mul(beta), p); // 9 DoubleMatrix temp = r; // 10 r = r2; r2 = temp;
Outlook
Better integration of different matrix classes.
Better performance through caching.
Support for sparse matrices and more compact storage schemes.
