August 16, 2012

I am fascinated by matrices. It is fundamental to so many computational algorithms these days, that I thought of having some fun writing a small JavaScript library about it. The source code is available for download from:

git clone https://github.com/gyaikhom/matrices.git

## Creates a matrix object with the supplied number of rows and columns.

Integer nrow Number of rows of the matrix.
Integer ncol Number of columns of the matrix.
Returns a matrix object of type Matrix
Matrix = function(nrow, ncol) {
var i, me = this;
me.r = nrow;
me.c = ncol;
me.v = new Array(nrow);
for (i = 0; i < nrow; ++i)
me.v[i] = new Array(ncol);
return me;
}

## Fills this matrix with values calculated using the supplied function.

Function f A function that accepts row and column indices to produce a value.
Returns the modified matrix of type Matrix
Matrix.prototype.fill = function(f) {
var i, j, me = this;
for (i = 0; i < me.r; ++i)
for (j = 0; j < me.c; ++j)
me.v[i][j] = f(i, j);
return me;
}

## Swaps two columns in this matrix.

By default, the matrix is modified after the swapping; however, this behaviour can be disabled.
Integer i The first column to swap.
Integer j The second column to swap.
Boolean m Should we leave the original matrix unmodified?
Returns a matrix with the columns swapped of type Matrix
Matrix.prototype.swapColumns = function(i, j, m) {
var me = this, r = undefined, k, t;
if (0 <= i && i < me.c && 0 <= j && j < me.c) {
r = m ? me.clone() : me;
for (k = 0; k < me.r; ++k) {
t = r.v[k][i];
r.v[k][i] = r.v[k][j];
r.v[k][j] = t;
}
}
return r;
}

## Swaps two rows in this matrix.

By default, the matrix is modified after the swapping; however, this behaviour can be disabled.
Integer i The first column to swap.
Integer j The second column to swap.
Boolean m Should we leave the original matrix unmodified?
Returns a matrix with the rows swapped of type Matrix
Matrix.prototype.swapRows = function(i, j, m) {
var me = this, r = undefined, k, t;
if (0 <= i && i < me.r && 0 <= j && j < me.r) {
r = m ? me.clone() : me;
for (k = 0; k < me.c; ++k) {
t = r.v[i][k];
r.v[i][k] = r.v[j][k];
r.v[j][k] = t;
}
}
return r;
}

## Clones this matrix.

Returns a clone of this matrix.
Matrix.prototype.clone = function() {
var me = this, r = new Matrix(me.r, me.c), i, j;
for (i = 0; i < me.r; ++i)
for (j = 0; j < me.c; ++j)
r.v[i][j] = me.v[i][j];
return r;
}

## Checks if the matrix is a square matrix.

Returns true if matrix is square; otherwise, false.
Matrix.prototype.isSquare = function() {
return this.c == this.r;
}

## Checks if this matrix equals the supplied matrix.

Matrix m A matrix to compare.
Returns true if the supplied matrix is equal to this matrix; otherwise, returns false.
Matrix.prototype.equals = function(m) {
var i, j, r = true, me = this;
if (m === null || !(m instanceof Matrix))
r = undefined;
else {
if (me.r == m.r && me.c == m.c)
for (i = 0; i < me.r; ++i) {
for (j = 0; j < me.c; ++j)
if (me.v[i][j] !== m.v[i][j]) {
r = false;
break;
}
if (!r) break;
}
else r = false;
}
return r;
}

## Checks if this matrix has the same order as the supplied matrix.

Matrix m A matrix to compare.
Returns undefined if the supplied object is not a matrix. Returns true if the supplied matrix has the same order as this matrix; otherwise, returns false.
Matrix.prototype.isSameOrder = function(m) {
var me = this;
return (m === null || !(m instanceof Matrix)) ? undefined :
(me.r == m.r && me.c == m.c) ? true : false;
}

## Prints this matrix to the browser console.

String prefix A prefix before the output.
String suffix A suffix after the output.
Returns the unmodified matrix.
Matrix.prototype.print = function(prefix, suffix) {
var i, j, me = this;
if (console) {
if (prefix) console.log(prefix);
for (i = 0; i < me.r; ++i)
for (j = 0; j < me.c; ++j)
console.log(me.v[i][j]);
if (suffix) console.log(suffix);
} else alert('Could not print matrix: console is undefined!');
return me;
}

## Returns the transpose of this matrix.

The original matrix is left unmodified.
Returns a new matrix which is the transpose.
Matrix.prototype.transpose = function() {
var i, j, me = this, t = new Matrix(me.c, me.r);
for (i = 0; i < me.r; ++i)
for (j = 0; j < me.c; ++j)
t.v[j][i] = me.v[i][j];
return t;
}

## Returns a scaled matrix after multiplying this matrix with the supplied scalar quantity.

By default, the original matrix will be modified; however, this behaviour can be disabled.
Real s The scalar quantity to multiply with.
Boolean m Should we leave the original matrix unmodified?
Returns a scaled matrix.
Matrix.prototype.scale = function(s, m) {
var i, j, r = m && m == true ? this.clone() : this;
for (i = 0; i < r.r; ++i)
for (j = 0; j < r.c; ++j)
r.v[i][j] *= s;
return r;
}

## Multiplies this matrix to the supplied matrix.

Both matrices are left unmodified. If this matrix is A and the supplied matrix is B and order is 'left', then the returned matrix is B x A; otherwise, this method returns A x B. By default, if the order is left unspecified, multiplication to the right, i.e., A x B, is assumed. This matrix and the supplied matrix is never modified.
Matrix m A matrix to multiply to.
String o Specifies the order of the multiplication.
Returns a matrix which is the product. If the supplied object is not a matrix, undefined is returned.
Matrix.prototype.multiply = function(m, o) {
var i, j, k, r = undefined, me = this, a, b;
if (m !== null && m instanceof Matrix) {
if (o && o == 'left') {
a = m;
b = me;
} else {
a = me;
b = m;
}
if (a.c == b.r) {
r = new Matrix(a.r, b.c);
for (i = 0; i < a.r; ++i)
for (j = 0; j < b.c; ++j) {
r.v[i][j] = 0;
for (k = 0; k < b.r; ++k)
r.v[i][j] += a.v[i][k] * b.v[k][j];
}
}
}
return r;
}

## Calculates the Hadamard product of this matrix and the supplied matrix.

Both matrices are left unmodified.
Matrix m A matrix to multiply with.
Returns a matrix which is the Hadamard product. If the supplied object is not a matrix, undefined is returned.
Matrix.prototype.hadamard = function(m) {
var i, j, r = undefined, me = this;
if (m !== null && m instanceof Matrix && me.isSameOrder(m)) {
r = new Matrix(m.r, m.c);
for (i = 0; i < m.r; ++i)
for (j = 0; j < m.c; ++j)
r.v[i][j] = m.v[i][j] * me.v[i][j];
}
return r;
}

## Calculates the Kronecker product of this matrix and the supplied matrix.

Both matrices are left unmodified. If this matrix is A and the supplied matrix is B and order is 'left', then the returned matrix is B x A; otherwise, this method returns A x B. By default, if the order is left unspecified, multiplication to the right, i.e., A x B, is assumed.
Matrix m A matrix to multiply to.
String o Specifies the order of the multiplication.
Returns a matrix which is the Kronecker product. If the supplied object is not a matrix, undefined is returned.
Matrix.prototype.kronecker = function(m, o) {
var i, j, k, l, r = undefined, me = this, a, b, x, y;
if (m !== null && m instanceof Matrix) {
r = new Matrix(me.r * m.r, me.c * m.c);
if (o && o == 'left') {
a = m;
b = me;
} else {
a = me;
b = m;
}
for (x = i = 0; i < a.r; ++i, x += b.r)
for (y = j = 0; j < a.c; ++j, y += b.c)
for (k = 0; k < b.r; ++k)
for (l = 0; l < b.c; ++l)
r.v[x + k][y + l] = a.v[i][j] * b.v[k][l];
}
return r;
}

## Returns a matrix which is the sum of this matrix and the supplied matrix.

By default, the original matrix will be modified; however, this behaviour can be disabled. The supplied matrix is never modified.
Matrix a A matrix to add.
Boolean m Should we leave the original matrix unmodified?
Returns a matrix which is the sum. If the supplied object is not a matrix, undefined is returned.
Matrix.prototype.add = function(a, m) {
var i, j, r = undefined, me = this;
if (me.isSameOrder(a)) {
r = m && m == true ? this.clone() : me;
for (i = 0; i < r.r; ++i)
for (j = 0; j < r.c; ++j)
r.v[i][j] += a.v[i][j];
}
return r;
}

## Returns a matrix which is the difference of this matrix and the supplied matrix.

By default, the original matrix will be modified; however, this behaviour can be disabled. The supplied matrix is never modified.
Matrix a A matrix to subtract from current matrix.
Boolean m Should we leave the original matrix unmodified?
Returns a matrix which is the difference. If the supplied object is not a matrix, undefined is returned.
Matrix.prototype.subtract = function(a, m) {
var i, j, r = undefined, me = this;
if (me.isSameOrder(a)) {
r = m && m == true ? this.clone() : me;
for (i = 0; i < r.r; ++i)
for (j = 0; j < r.c; ++j)
r.v[i][j] -= a.v[i][j];
}
return r;
}

## Creates an identity matrix object with the supplied order.

Integer o Order of the matrix.
Returns a matrix object.
Matrix.prototype.Identity = function(o) {
var i, j, m = new Matrix(o, o);
for (i = 0; i < o; ++i) {
for (j = 0; j < o; ++j)
m.v[i][j] = 0;
m.v[i][i] = 1;
}
return m;
}

## Creates an square matrix object with the supplied order.

Integer o Order of the matrix.
Returns a matrix object.
Matrix.prototype.Square = function(o) {
var i, j, m = new Matrix(o, o);
for (i = 0; i < o; ++i)
for (j = 0; j < o; ++j)
m.v[i][j] = 0;
return m;
}

## Calculates the trace of this matrix, if it is square.

Returns the trace of the square matrix, or undefined, if not square.
Matrix.prototype.trace = function() {
var r = undefined, i, me = this;
if (me.isSquare())
for (r = 0, i = 0; i < me.r; ++i)
r += me.v[i][i];
return r;
}

## Checks if the matrix is a null matrix.

Returns true if the matrix is null; otherwise, false is returned.
Matrix.prototype.isNull = function() {
var i, j, me = this, r = true;
for (i = 0; i < me.r; ++i) {
for (j = 0; j < me.c; ++j)
if (me.v[i][j] != 0) {
r = false;
break;
}
if (!r) break;
}
return r;
}

## Checks if the matrix is a zero matrix (which is the same as being null).

Returns true if the matrix is a zero matrix; otherwise, false is returned.
Matrix.prototype.isZero = Matrix.prototype.isNull;

## Checks if the matrix is a diagonal matrix.

Returns true if the matrix is diagonal; otherwise, false is returned.
Matrix.prototype.isDiagonal = function() {
var i, j, me = this, r = false;
if (me.isSquare()) {
r = true;
for (i = 0; i < me.r; ++i) {
for (j = 0; j < me.c; ++j)
if (i != j && me.v[i][j] != 0) {
r = false;
break;
}
if (!r) break;
}
}
return r;
}

## Checks if this matrix is upper triangular.

Returns if square and upper triangular, return true; otherwise, return false.
Matrix.prototype.isUpperTriangular = function() {
var r = true, i, j, me = this;
if (me.isSquare()) {
for (i = 0; i < me.r; ++i) {
for (j = 0; j < i; ++j)
if (me.v[i][j] !== 0) {
r = false;
break;
}
if (!r) break;
}
} else r = false;
return r;
}

## Checks if this matrix is lower triangular.

Returns if square and lower triangular, return true; otherwise, return false.
Matrix.prototype.isLowerTriangular = function() {
var r = true, i, j, me = this;
if (me.isSquare()) {
for (i = 0; i < me.r; ++i) {
for (j = i + 1; j < me.c; ++j)
if (me.v[i][j] !== 0) {
r = false;
break;
}
if (!r) break;
}
} else r = false;
return r;
}

## Checks if this matrix is symmetric.

Returns if square and symmetric, return true; otherwise, return false.
Matrix.prototype.isSymmetric = function() {
var r = true, i, j, me = this;
if (me.isSquare()) {
for (i = 0; i < me.r; ++i) {
for (j = 0; j < me.c; ++j)
if (me.v[i][j] !== me.v[j][i]) {
r = false;
break;
}
if (!r) break;
}
} else r = false;
return r;
}

## Checks if this matrix is skew-symmetric.

Returns if square and skew-symmetric, return true; otherwise, return false.
Matrix.prototype.isSkewSymmetric = function() {
var r = true, i, j, me = this;
if (me.isSquare()) {
for (i = 0; i < me.r; ++i) {
for (j = 0; j < me.c; ++j)
if (me.v[i][j] != -me.v[j][i]) {
r = false;
break;
}
if (!r) break;
}
} else r = false;
return r;
}

## Extracts lower/upper triangular factors from this matrix, which stores a packed LU decomposition result (L as unit lower triangular matrix).

Boolean l To extract the lower triangular matrix, pass true; otherwise, pass false.
Returns a lower/upper triangular matrix. If the LU decomposition is not a matrix, undefined is returned.
Matrix.prototype.luExtract = function(l) {
var i, j, me = this, order = me.r, r = me.clone(), v = r.v;
if (l)
for (i = 0; i < order; ++i) {
for (j = i + 1; j < order; ++j)
v[i][j] = 0;
v[i][i] = 1;
}
else
for (i = 0; i < order; ++i)
for (j = 0; j < i; ++j)
v[i][j] = 0;
return r;
}

## Separates the LU decomposition result into L and U factors.

Return an object that contains two properties L and U, where, each respectively gives the L and U decomposition of this matrix is returned. If the matrix is not square, or LU decomposition not stable, undefined is returned.
Matrix.prototype.luSeparate = function() {
var me = this;
return {
L: me.luExtract(true), /* extract lower triangular factor */
U: me.luExtract(false) /* extract upper triangular factor */
};
}

## Calculates the LU decomposition of this matrix without pivoting.

We use in-place decomposition, and the factors are separated out into separate L and U matrices before they are returned.

This separation into two matrices can be disabled, in which case, both L and U are placed in the same matrix, with the understanding that the diagonal of L is all 1s.

This method leaves the current matrix unmodified.

Boolean s Should we disable separation of the L and U matrices? By default, do separate.
If separation is disabled, a single matrix with both L and U are returned. If separation is enabled, an object that contains two properties L and U, where, each respectively gives the L and U decomposition of this matrix is returned. If the matrix is not square, or LU decomposition not stable, undefined is returned.
Matrix.prototype.lu = function(s) {
var me = this, r = undefined;
if (me.isSquare()) {
var t = me.clone(), i, j, k, x = me.r, y = x - 1, v = t.v;
for (k = 0; k < y; ++k) {
/* column normalisation */
for (i = k + 1; i < x; ++i)
v[i][k] = v[i][k] / v[k][k];

/* submatrix modification */
for (i = k + 1; i < x; ++i)
for (j = k + 1; j < x; ++j)
v[i][j] = v[i][j] - v[i][k] * v[k][j];
}
r = s && s == true ? t : t.luSeparate();
}
return r;
}

## Calculates the LU decomposition of this matrix using Doolittle algorithm.

We use in-place decomposition, and the factors are separated out into separate L and U matrices before they are returned.

This separation into two matrices can be disabled, in which case, both L and U are placed in the same matrix, with the understanding that the diagonal of L is all 1s.

This method leaves the current matrix unmodified.

Boolean s Should we disable separation of the L and U matrices? By default, do separate.
If separation is disabled, a single matrix with both L and U are returned. If separation is enabled, an object that contains two properties L and U, where, each respectively gives the L and U decomposition of this matrix is returned. If the matrix is not square, or LU decomposition not stable, undefined is returned.
Matrix.prototype.doolittle = function(s) {
var me = this, r = undefined;
if (me.isSquare()) {
var t = me.clone(), i, j, k, order = me.r, v = t.v;
for (i = 0; i < order; ++i) {
/* determine upper triangular factor */
for (j = i; j < order; ++j)
for (k = 0; k < i; ++k)
v[i][j] -= v[i][k] * v[k][j];
if (v[i][i] == 0) {
t = undefined; /* this is a singular matrix */
break;
}
/* determine lower triangular factor */
for (j = i + 1; j < order; ++j) {
for (k = 0; k < i; ++k)
v[j][i] -= v[j][k] * v[k][i];
v[j][i] /= v[i][i];
}
}
r = s && s == true ? t : t.luSeparate();
}
return r;
}

## Calculates the Cholesky decomposition of this matrix.

If the matrix A is symmetric and positive-definite, this method returns a matrix L such that A = LL*, where L is a lower triangular matrix with a strictly positive diagonal, and L is the transpose of L. This method does not modify the matrix.

For now, we only deal with real numbers, and not imaginary numbers, in which case L* should be the conjugate transpose of L.

If symmetric and positive-definite, L is returned; otherwise, undefined is returned.
Matrix.prototype.cholesky = function() {
var me = this, r = undefined;
if (me.isSymmetric()) {
var i, j, k, x, y, s, order = me.r, v = me.v, w;
r = new g7.squarematrix(order);
w = r.v;
for (i = 0; r && i < order; ++i) {
x = i + 1;
for (j = 0; r && j < x; ++j) {
y = 0;
for (k = 0; k < j; ++k)
y += w[i][k] * w[j][k];
if (i == j) {
s = v[i][i] - y;
if (s < 0) {
throw new Error('Matrix is not positive-definite');
r = undefined;
} else w[i][i] = Math.sqrt(s);
} else w[i][j] = 1 / w[j][j] * (v[i][j] - y);
}
}
}
return r;
}