Something about matrices

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

Table of contents

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;
}
comments powered by Disqus