vkostyukov / la4j Goto Github PK

Shouldn't

private static class DecMatrixFunction
implements MatrixFunction {
@OverRide
public double evaluate(int i, int j, double value) {
return value + 1.0;
}
}

Read

private static class DecMatrixFunction
implements MatrixFunction {
@OverRide
public double evaluate(int i, int j, double value) {
return value - 1.0;
}
}

instead? You don't seem to be using it anywhere at the moment

This engancement was requested by Xiaorui Jiang. Here is his reference implementation:

Matrix[][] createBlockMatrix(A, B, C, D) {
 Matrix[][] blocks = new Matrix[2][2];
 blocks[0][0] = A;
 blocks[0][1] = B;
 blocks[1][0] = C; 
 blocks[1][1] = D;
 return blocks;
}

So, createBlockMatrix method should be vararg like this:

Matrix[][] crateBlockMatrix(Matrix ... matrices);

The update() method has default implementation in AbstractMatrix class. Actually, it is enough for dense matrices but we can reduce ovehead by implementing it for sparse matices.

Here is the reference implementation https://sites.google.com/site/indy256/algo/determinant

Something like:

Matrix resizeColumns(int i);
Matrix resizeRows(int i);

See subj.

There are two major test classes AbstractMatrixTest and AbstractVectorTess. So, we have to refactor these classes to use different input data. The main problem is these classes use only 3x3 matrices, but somtimes we need to check other size.

There is prove of concept for new teset methods:

testTranspose_3x3 - for methods w/o paramters | transposes 3x3 matrix;

testMultiply_3x3_3x3 - for methods w/parameters | muptiplies two matrices with size 3x3;

testMultiply_3x3_double - for methods w/parameters | muptiplies 3x3 matrix to double;

testMultiply_3x3_3 - for methods w/parameters | muptiplies 3x3 matrix to 3-items vector;

and so on.

We need better test coverage for rectagle (n>m, n<m) and square matrices.

Here is draft: http://arxiv.org/pdf/1209.3995.pdf

This is an inovative algorithgm designed in 2012 by Joerg Fliege. It would be awesome if la4j provided a support for this approach.

BTW, Gaussian Ellimination (that uses by default in la4j) works for O(n^3)

The class SymmetricMatrixPredicate should be improved. There is a mamory and performance issues.

It could be a 0 (zero) as parameter for matrix.set(). So, we have to remove the existed cell in this case.

Here is the right place to put the implementation

for (int jj = columnPointers[j]; jj < columnPointers[j + 1]; jj++) {
            if (rowIndices[jj] == i) {

                // TODO: if value is zero:
                // clear the value cell

                values[jj] = value;
                return;
            }
        }

As both operands are integers, an integer division is performed:

public double density() {
    return cardinality / length;
}

Thus, CompressedVector.density() always returns 0.0.

Update: Same issue with AbstractCompressedMatrix.density().

Other words, the resize() method should return new matrix.

There are a lot of articles about linear algebra at WlframAlpha:http://mathworld.wolfram.com/topics/Matrices.html. It would be great if each la4j class had a reference to WA article with simple note aout parameters and exceptions.

Matrices in la4j use 0-based indexing of rows and columns while Matrix Market format requires matrix indices to be 1-based (Matlab like).

Matrix Market Format Specifications - http://math.nist.gov/MatrixMarket/formats.html#MMformat

p.s.: thank you for your great work!

There is each() method that should iterate through the all elements of entry. It would be great to have separate method through iterating through non-zero values. Like this:

The Matrix should contans following methods:

void each(MatrixProcedure procedure); //
void eachInRow(int i, MatrixProcudure procedure); // iterating row elements
void eachInColumn(int j, MatrixProcudure procedure); // iterating column elements

void eachnz(MatrixProcedure procedure); // iterating through non-zero elements
void eachnzInRow(int i, MatrixProcedure procedure); // 
void eachnzInColumn(int i, MatrixProcedure procedure); // 

The Vector:

void each(VectorProcedure procedure); //
void eachnz(VectorProcedure procedure); //

Class name IdentityMattixSource contains a typing error.

Here is the mentioned operator:

return new Matrix[] { q.multiply(-1), r };

We have to change the algorithm to remove this operation.

Should be something like following:

Matrix a = new CRSMatrix(...)
Matrix b = a.shuffle();

So it has the signature in Matrix class:

Matrix shuffle();

There are sparse entries in la4j, like CRSMatrix and CCSMatrix. So, it is not necessary to iterate through the all elements of these matrices while multiplying (for example). So, it would be great to develop proof-of-concept that will allow to hide details from developes and using the same API provide a more efficient way to perform basic operations.

In terms of matrix multiplying, it should be:

Matrix a = new CRSMatrix(...);
Matrix b = new Basic2DMatrix(...);

Matrix c = a.multiply(b);

Just an Idea: Probably, it would be the good idea to incapsulate all operations in special objects, which know the internal structure (or some other details) and generate such objects in matrices methods like multiply from two matrix operands.

I'm open for descussion.

Subj. Also we need to think about using these formats with dense matrices.

double sum();

double product();

And the example:

Matrix a = new CRSMatrix(...);

double d1 = a.sum();
double d2 = a.product();

I've found this implementation in numbers.js library. Here is the code:

if (arrA[0].length === arrB.length) {
  var result = new Array(arrA.length);

for (var x = 0; x < arrA.length; x++) {
  result[x] = new Array(arrB[0].length);
}

var arrB_T = transpose(arrB);

for (var i = 0; i < result.length; i++) {
  for (var j = 0; j < result[i].length; j++) {
    result[i][j] = dotproduct(arrA[i],arrB_T[j]);
  }
}
return result;

} else {
throw new Error("Array mismatch");
}

Probabbly, it works faster than current implementation.

There are code blocks like following in la4j:

a.set(i, j, a.get(i, j) * 10);

It can be replaced with update() method:

a.update(i, j, Matrices.asMulFunction(10));

The debug matrices should log out information about internal proccess like multiplication. For example toString() method of debug matrix should out:

CRS = { values = [1, 2, 3]; columnIndices = [4, 5, 6]; rowPointers = [7, 8, 9] }

Probabbly, it can be used instead of compound operators like following. It would be nice to have the similar transform methods:

Matrix transform(int i, int j, MatrixFunction function);

void Matrix.update(int i, int j, MatrixFunction function);
void Matrix.update(MatrixFunction function);

Matrices.asPlusFunction(10); // x = x + y
Matrices.asMinusFunction(10); // x = x - y
Matrices.asMulFunction(10); // x = x * y
Matrices.asDivFunction(10); // x = x / y

And the same for vectors.

The difference between transform and update is that the first one is returning new matrix but the second one performing changes In-Place.

So, it shoud add some overhead but I belive in JIT and hope that it will be inlined ad will cost 0.

This issue is submitted by Daniel Renshaw:

There are two type of products which should be calculated for vectros:

double innerProduct(Vector vector); //multiplies a row-vector by column-vector and produces a single value
Matrix outerProduct(Vector vector); //multiplies a column-vector by row-vector and produces a matrix

See here for details: http://lpsa.swarthmore.edu/BackGround/RevMat/MatrixReview.html

I'm unable to eigendecompose some sparse 24x24 symmetric matrices. The following example demonstrates the problem using la4j 0.4.0 (the most current version).

import org.la4j.matrix.Matrices;
import org.la4j.matrix.Matrix;

public class Program {
    public static void main(String[] args) {
        Matrix matrix = Matrices.BASIC2D_FACTORY.createSquareMatrix(24);
        matrix.set(0, 16, 6.488917882340233);
        matrix.set(1, 14, 1.0500849109563222);
        matrix.set(1, 18, 8.057338506021317);
        matrix.set(2, 18, 5.600914629504165);
        matrix.set(4, 22, 9.912502191842943);
        matrix.set(5, 13, 10.983197714528727);
        matrix.set(8, 13, 6.5109096798879476);
        matrix.set(9, 16, 1.321253910633423);
        matrix.set(10, 14, 7.175174747019432);
        matrix.set(10, 17, 3.0219813644087403);
        matrix.set(11, 18, 3.5278266874188975);
        matrix.set(12, 13, 7.425806151241948);
        matrix.set(12, 15, 3.7347574055072483);
        matrix.set(13, 5, 10.983197714528727);
        matrix.set(13, 8, 6.5109096798879476);
        matrix.set(13, 12, 7.425806151241948);
        matrix.set(14, 1, 1.0500849109563222);
        matrix.set(14, 10, 7.175174747019432);
        matrix.set(15, 12, 3.7347574055072483);
        matrix.set(15, 22, 3.1315457717146056);
        matrix.set(16, 0, 6.488917882340233);
        matrix.set(16, 9, 1.321253910633423);
        matrix.set(17, 10, 3.0219813644087403);
        matrix.set(17, 18, 8.059047125400168);
        matrix.set(18, 1, 8.057338506021317);
        matrix.set(18, 2, 5.600914629504165);
        matrix.set(18, 11, 3.5278266874188975);
        matrix.set(18, 17, 8.059047125400168);
        matrix.set(22, 4, 9.912502191842943);
        matrix.set(22, 15, 3.1315457717146056);
        System.out.println(matrix);

        // The following line does not complete successfully.
        Matrix[] vd = matrix.decompose(Matrices.EIGEN_DECOMPOSITOR);

        System.out.println(vd[0]);
        System.out.println(vd[1]);
    }
}

This matrix can be decomposed by Matlab fine. To see this, save the following as eigtestmat.tsv and execute the following Matlab code.

24  24  0
1   17  6.48891788234023
2   15  1.05008491095632
2   19  8.05733850602132
3   19  5.60091462950416
5   23  9.91250219184294
6   14  10.9831977145287
9   14  6.51090967988795
10  17  1.32125391063342
11  15  7.17517474701943
11  18  3.02198136440874
12  19  3.5278266874189
13  14  7.42580615124195
13  16  3.73475740550725
14  6   10.9831977145287
14  9   6.51090967988795
14  13  7.42580615124195
15  2   1.05008491095632
15  11  7.17517474701943
16  13  3.73475740550725
16  23  3.13154577171461
17  1   6.48891788234023
17  10  1.32125391063342
18  11  3.02198136440874
18  19  8.05904712540017
19  2   8.05733850602132
19  3   5.60091462950416
19  12  3.5278266874189
19  18  8.05904712540017
23  5   9.91250219184294
23  16  3.13154577171461

load eigtestmat.tsv
X=spconvert(eigtestmat)
[V,D]=eigs(X)

It's not just this matrix; I've tried another, of similar structure, which exhibits the same problem.

Ths method (in pair of transform()) should implement "map-reduce" paradigm. So, matrices aready has 'map' equivalent such transrom(). The fold() should be reduce() equivalent.

For matrices:

double fold(MatrixFunction function);

For vectors:

double fold(VectorFunction function);

Subj.

Matrix createConstantMatrix(rows, columns, value);

Reference implemenation in Scala

def rotateOutOfPlace(a: Array[Array[Int]]): Array[Array[Int]]= {
  val b: Array[Array[Int]] = Array.fill(a(0).length) { new Array(a.length) }

  for (i <- 0 until a.length) {
    for (j <- 0 until a(i).length) {
      b(j)(a.length - i -1) = a(i)(j) 
    }
  }

  return b
}

assert { rotateOutOfPlace(Array(Array(1 ,2), Array(3, 4))).deep == Array(Array(3, 1), Array(4, 2)).deep }

About positive defenite matrix: http://mathworld.wolfram.com/PositiveDefiniteMatrix.html

Reference implementation:

 private boolean isPositiveDefinite(Matrix matrix) {

        //
        // TODO: create a MatrixPredicate for it
        //

        int n = matrix.rows();
        boolean result = true;

        Matrix l = matrix.blank();

        for (int j = 0; j < n; j++) {

            Vector rowj = l.getRow(j);

            double d = 0.0;
            for (int k = 0; k < j; k++) {

                Vector rowk = l.getRow(k);
                double s = 0.0;

                for (int i = 0; i < k; i++) {
                    s += rowk.unsafe_get(i) * rowj.unsafe_get(i);
                }

                s = (matrix.unsafe_get(j, k) - s) / l.unsafe_get(k, k);

                rowj.unsafe_set(k, s);
                l.setRow(j, rowj);

                d = d + s * s;
            }

            d = matrix.unsafe_get(j, j) - d;

            result = result && (d > 0.0);

            l.unsafe_set(j, j, Math.sqrt(Math.max(d, 0.0)));

            for (int k = j + 1; k < n; k++) {
                l.unsafe_set(j, k, 0.0);
            }
        }

        return result;
    }

Here is the defenition http://en.wikipedia.org/wiki/Kronecker_product

Here is the contract for new method:

Matrix kronecker(Matrix matrix);

Hi,
Could you tell me if there is some documentation (like Javadoc) other than the examples & samples in the project description ?
Thanks !

CompressedVector.swap(int,int) is broken.

new CompressedVector(new double[]{0.0, 1.0, 0.0, 2.0}).swap(1, 2);

Expected result:

[0.0, 0.0, 1.0, 2.0]

Actual result:

[0.0, 0.0, 1.0, 0.0]

There seems to be a bug in la4j's eigenvalue calculations. I ran the following code:

public static void main(String args[])
{
    Factory crs = new CRSFactory();
    Matrix matrix = crs.createRandomSymmetricMatrix(6);
    Matrix[] eig = matrix.decompose(Matrices.EIGEN_DECOMPOSITOR,Matrices.CRS_FACTORY);
}

About half the time I run this code it terminates quickly (and gives the right output if I add a few lines to display it), but about half the time it seems to hang up in an infinite loop. For larger matrices it almost always hangs up.

Matrices should have an ability to calc the rank. Here is rank's explanation: http://mathworld.wolfram.com/MatrixRank.html

la4j can calculate determinant for any matrix. It is wrong, since only square matrices has determinant property.

It will allow to remove ugly unsafe_ methods. But I have to think about backward capability. For instance: current matrix classes should be safe.

Read the subj.

Here is some proof-of-concept for new methods.

In Matrix:

void add(int i, int j, double value);
void multiply(int i, int j, double value);

And the same ones for Vector.

Here is the possible contract:

Matrix.java

Matrix slice(int i, int j, int k, int l);
Matrix sliceLeft(int i, int j);
Matrix sliceRight(int i, int j);

Vector.java

Vector slice(int i, int j);
Matrix sliceLeft(int i);
Matrix sliceRight(int i);

http://en.wikipedia.org/wiki/Hadamard_product_(matrices)

Here is the contract

Matrix hadamardProduct(Matrix matrix);

Matrix hadamardProduct(Matrix matrix, Factory factory);

See subj.

double diagonalProduct();

I think that function .getArray() should be implemented both for Vector and Matrix.

AbstractVector.subtract() throws a NullPointerException, if its first argument is null.

See

There is MAX_ITEARATIONS in Eigen

public static final int MAX_ITERATIONS = 10000000;

We need to make sure that this values is big enough. Probably, it might be interesting to compile some table with experemental values number of performed iterations to decompose the matrix with size MxN. For diferent matrix sizes (small and big).

Current implementation allows to get i-th element of compressed vector in linear time O(n) where n is size of vector. There is a low-level-fruit of improving this time. The idea is to use ordered storage and binary search in get and set operations. It will allow to reduce get/set operations time from linear to logarithmic O(log n)

There are several methods that generates new matrices (random matrix, indentity, symmetric). So, these implmenations should use the "raw" matrix structure to improve the performance of this code.