Friday, March 6, 2015

Matrices calculations for numerical analysis and linear algebra in Java.

Matrices.java
This is a class containing implementations of known numerical analysis and linear algebra algorithms about matrices.Check github for examples.

Algorithms:

Google's PageRank implementation in Java

Calculating the Google matrix from an adjacency matrix:

public double[][] getGoogleMatrix(int[][] a,double q) {

    int n=a.length;
    double[][] g = new double[n][n];
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            g[i][j] = ((q / n) + ((a[j][i] * (1 - q)) / sumOfLine(a[j])));
        }
    }
    return g;
}

Ranking the nodes can be done using the values from the normalized eigenvector of the google matrix:
public double[] normalize(double[] a){

    double s=0;
    for (int i = 0; i < a.length; i++) {
        s+=abs(a[i]);
    }
    double[] n;
    n=multiply(a,1/s);
    for (int i = 0; i < n.length; i++){
        n[i]=(double) round(n[i] * 100000000) / 100000000;
    }
    return n;
}

Power Iteration method to calculate the eigenvector.

This is part of the Matrices project.Refer to Basic matrix calculations in Java for the missing secondary methods.

Power Iteration for the largest eigenvector in Java

Wikipedia

public double[] powerIt(double[][] a) {

    int n=a.length;
    double e = 0.5 * pow(10, -7);
    double[] b = new double[n];
    Random r = new Random(System.currentTimeMillis());
    for (int i = 0; i < n; i++) {
        b[i] = r.nextDouble();
    }

    double l0;
    double l1 = firstNonZero(b);
    do {
        l0 = l1;
        b = multiply(a, b);
        l1 = firstNonZero(b);
        b=multiply(b,1/l1);
    } while (abs(l1 - l0) > e);
    return b;
}

This is part of the Matrices project.Refer to Basic matrix calculations in Java for the missing secondary methods.

Gauss Seidel method in Java

Wikipedia

double s1;

double s2;
int loop;
double[] x;
public double[] gaussSeidel(final double[][] a,double[] b,double[] start){
    double e=0.5*pow(10,-6);
    x=copy(start);
    final int n=start.length;
    double[] prev;
    do{
        prev=copy(x);
        for (loop = 0; loop < n; loop++) {
            s1=0;
            s2=0;
            Thread t1 = new Thread() {
                @Override
                public void run() {
                    for (int j = 0; j < loop; j++) {
                        s1 += a[loop][j] * x[j];
                    }
                }
            };
            Thread t2 = new Thread() {
                @Override
                public void run() {
                    for (int k = loop + 1; k < n; k++) {
                        s2 += a[loop][k] * x[k];
                    }
                }
            };
            t1.start();
            t2.start();
            try {
                t1.join();
                t2.join();
            } catch (InterruptedException ex) {

            }
            x[loop] = (1.0 / a[loop][loop]) * (b[loop]-s1-s2);
        }
    }while(normInf(dif(x,prev))>=e);
    return x;
}

This is part of the Matrices project.Refer to Basic matrix calculations in Java for the missing secondary methods.

Cholesky decomposition in Java

Wikipedia

public double[][] cholesky(double[][] a){

    double[][] l=new double[a.length][a.length];
    double s;
    for (int i = 0; i < a.length; i++) {
        for (int j = 0; j < i; j++) {
            s=0;
            for (int k = 0; k < j; k++) {
                s+=l[i][k]*l[j][k];
            }
            l[i][j]=(a[i][j]-s)/l[j][j];
        }
        s=0;
        for (int k = 0; k < i; k++) {
            s+= pow(l[i][k],2);
        }
        l[i][i]=sqrt((a[i][i]-s));
    }
    return l;
}

This is part of the Matrices project.

Gaussian ellimination using LU factorization in Java

public double[] gauss(double[][] a,double b[]){

    int n=a.length;
    double[][] p=identity(n);
    double[][] l=identity(n);
    double[][] u=pivot(a,p);
    for (int i = 1; i < n; i++) {
        for (int j = 0; j < i; j++) {
            if(u[i][j]!=0){
                l[i][j]=u[i][j]/u[j][j];
                u[i]=multiplyAndSumLines(u[j],u[i],-l[i][j]);
            }
        }
    }
    double[] y=new double[n];
    double[] z=multiply(p,b);
    double s;
    for (int i = 0; i < n; i++) {
        s=0;
        for (int j = 0; j < i; j++) {
            s+=l[i][j]*y[j];
        }
        y[i]=(z[i]-s)/l[i][i];
    }
    double[] x=new double[n];
    for (int i = n-1; i >= 0; i--) {
        s=0;
        for (int j = n-1; j > i; j--) {
            s+=u[i][j]*x[j];
        }
        x[i]=(y[i]-s)/u[i][i];
    }
    return x;
}

public double[][] pivot(double[][] a, double[][] p) {

    double[][] acopy=copy(a);
    for (int i = 0; i < a.length; i++) {
        int imax=jMaxOf(getColumn(acopy,i),i);
        if(imax!=i){
            double temp;
            for (int j = 0; j < a.length; j++) {
                temp=p[i][j];
                p[i][j]=p[imax][j];
                p[imax][j]=temp;

                temp=acopy[i][j];
                acopy[i][j]=acopy[imax][j];
                acopy[imax][j]=temp;
            }
        }
    }
    return acopy;
}

This is part of the Matrices project.Refer to Basic matrix calculations in Java for the missing secondary methods.

Basic matrix calculations in Java

MatricesBasic.java
A set of methods used in the Matrices project.

Commonly used:

public double[][] identity(int n){

    double[][] array=new double[n][n];
    for (int i = 0; i < n; i++) {
        array[i][i]=1;
    }
    return array;
}
public double[] multiply(double[] a,double d){
    double[] r=new double[a.length];
    for (int i = 0; i < a.length; i++) {
        r[i]=(a[i]*d);
    }
    return r;
}
public double[] multiplyAndSumLines(double[] a,double[] b,double m){
    double[] r=multiply(a,m);
    for (int i = 0; i < a.length; i++) {
        r[i]+=b[i];
    }
    return r;
}
public double normInf(double[] a){
    double max= abs(a[0]);
    for (int i = 1; i < a.length; i++) {
        if(abs(a[i])>max)
            max=abs(a[i]);
    }
    return max;
}
public int jMaxOf(double[] A,int from) {
    int mj=from;
    double m= abs(A[from]);
    for (int j = from; j < A.length; j++) {
        if(abs(A[j])>m){
            m= abs(A[j]);
            mj=j;
        }
    }
    return mj;
}