[FMX] Matrix algebra and matrix inverse

[chris23]

Has anyone here ever developed a function (or general library of relevant functions) that will allow one to find the inverse of a matrix? I've searched the forum archives, as well as the internet more generally, but haven't found anything. If anyone has some leads, I hope you'll pass them along.

If I don't hear anything within the next few days, I might try to develop something on my own, but I'd hate to do so if someone has already done the dirty work.

A solution to this problem would be very useful for 3-D work or anyone doing mathematical programming more generally.

Thanks,
Chris

[eyezberg]

I'm not into that at all, check here maybe:
http://www.layer51.com/proto/

[chris23]

Hi, eyezberg. I've already checked that site. They do have some nice functions for calculating the transpose of a matrix and other operations, but no functions designed for calculating the inverse of a matrix.

[eyezberg]

As I said, not into that stuff at all, you might want to ask Robert Penner or Guy Watson (Flashguru.co.uk)...

[chris23]

Ok, I finally figured out a way to solve for the inverse of a matrix using FMX ActionScript. It took a while, but I think it will be worth it. This will be very helpful in some of my applications, and I hope there may be others out there who will find it useful too.

The code is listed in a separate post below in order to deal with the Kirupa Forum message length limitations. I'll provide a brief overview of a problem that might motivate the need to find the inverse of a matrix or multidimensional array in this post. In the next post I'll provide some examples of how the matrix prototypes can be used. In the final post I'll post the codes.

Assume we were working with the following system of equations:

1[x1] + 2[x2] + 3[x3] = 14
1[x1] + 3[x2] + 2[x3] = 13
2[x1] + 4[x2] + 1[x3] = 13

This set of equations can be written in matrix form as:

   A     X   =   Y
[1 2 3] [x1]   [14]
[1 3 2] [x2]   [13]
[2 4 1] [x3]   [13],

where A is a 3 x 3 matrix of weights, X is a 3 x 1 matrix of unknowns, and Y is a 3 x 1 matrix containing outcomes.

According to this matrix equation (i.e., AX = Y), some set of outcomes, Y, can be modeled as a weighted function of X. The trick is to solve this system of simultaneous equations for X.

To do so it is necessary to find the inverse of matrix A. The inverse of a matrix behaves like a reciprocal in regular algebra. For example, if we were working with non-matrices and had the equation 2*x = 6, we could solve the equation for x by multiplying both sides of the equation by the reciprocal of 2 (i.e., .5). Doing so would yield .5(2)*x = (.5)6 -> 1*x = 3 -> x = 3. Calculating an inverse of a matrix, however, is fairly complex. Once it has been calculated, however, the same basic logic can be used to solve the equations. In the matrix example above:

AX = Y 
inverse(A)AX = inverse(A)Y
IX = inverse(A)Y
X = inverse(A)Y,

where I is an identity matrix--a matrix of zeros with 1's on the diagonal. (An identity matrix behaves like a 1 in regular algebra.)

In this example, the inverse of matrix A is as follows:

[  1 -2   1]
[-.6  1 -.2]
[ .4  0 -.2].

We can now use this inverse to solve the system of equations for X:

 inverse(A)     A     X   =   inverse(A)  Y
[  1 -2   1] [1 2 3] [x1]   [  1 -2   1] [14]
[-.6  1 -.2] [1 3 2] [x2]   [-.6  1 -.2] [13]
[ .4  0 -.2] [2 4 1] [x3]   [ .4  0 -.2] [13],

   I     X   =   inverse(A)  Y
[1 0 0] [x1]   [  1 -2   1] [14]
[0 1 0] [x2]   [-.6  1 -.2] [13]
[0 0 1] [x3]   [ .4  0 -.2] [13],

  X   = inverse(A)Y
[x1]    [1]
[x2]    [2]
[x3]    [3].

Substituting this solution into the equation suggests that

   A     X   =   Y
[1 2 3] [1]   [14]
[1 3 2] [2]   [13]
[2 4 1] [3]   [13].

[chris23]

I've pasted the code below for a function that will find the inverse of a matrix. In the process of developing this function, however, I found it useful to develop a broader array of functions that can be used to perform matrix operations. Instead of creating a new class of matrix objects, these functions simply allow matrix operations to be performed on multidimensional arrays. In other words, I've created several array prototypes that can be used to not only find the inverse of a multidimensional array, but perform other matrix functions as well (i.e., matrix multiplication, matrix transpose).

Since these are prototypes, they probably need very little explanation. Nonetheless, here is a set of commands that can be used to generate the numbers described above:

// Create a 3 x 3 matrix and insert each element one by one
// (note this is obviously not the most efficient way to enter in the values of
// a multidimensional array, but it is useful here for illustrating the various 
// functions
A = Matrix(3, 3);
A[0][0] = 1;
A[1][0] = 1;
A[2][0] = 2;
A[0][1] = 2;
A[1][1] = 3;
A[2][1] = 4;
A[0][2] = 3;
A[1][2] = 2;
A[2][2] = 1;
//
trace("Display matrix A");
A.print();
//
trace("Display matrix/vector Y")
Y = new Array(14, 13, 13);
//
trace("Find and Display Ai -- the inverse of A");
Ai = A.inverse();
Ai.print();
//
// Solve for X 
// X = (Ai)(Y)
trace("Solve the equations for X and display the solution ");
X = Ai.matrixVectorMult(Y);
X.print();

Here is the output of the trace commands:

Display matrix A
[1,2,3]
[1,3,2]
[2,4,1]

Display matrix/vector Y
Find and Display Ai -- the inverse of A
determinant = -5
[1,-2,1]
[-0.6,1,-0.2]
[0.4,0,-0.2]

Solve the equations for X and display the solution
[0.999999999999996]
[2]
[3]

The library of prototypes are included in the next post. Some of the functions are not terribly efficient. If anyone has suggestions on how to make all of this work more smoothly, I hope you'll post them. As I continue to work on this stuff, I'll post modifications to the code.

--Chris

[chris23]

Just copy and paste this code into the first frame of the movie. Additional description is provided in comments. I'll modify these functions within the post as I have time to clean them.
--Chris

//
// Matrix Algebra Library for Flash MX ActionScript
//   The functions can be used to apply matrix-based operations to arrays. 
//   No new objects or classes needed.
//
// R. Chris Fraley
// April 2003
//
// -------------------------------------------
// Creating and printing matrices
// -------------------------------------------
//
// Create a new array with matrix structure
//  (matrix will be null by default)
//
function Matrix(nrows, ncols) {
	var output = new Array();
	if (nrows != 0 && ncols != 0) {
		for (var i = 0; i < nrows; i++) {
			var temp2 = new Array();
			for (var j = 0; j < ncols; j++) {
				temp2.push(0);
			}
			output[i] = temp2;
		}
	}
	if (nrows == 1) {
		var temp2 = new Array();
		for (var j = 0; j < ncols; j++) {
			temp2.push(0);
		}
		output[0] = temp2;
	}
	if (ncols == 1) {
		var temp2 = new Array();
		for (var j = 0; j < nrows; j++) {
			//temp2.push(0);
			output[j] = Array();
			output[j] = 0;
		}
	}
	output.dim();
	return (output);
}
//
// Run this function to assign nrow and ncol properties to an array
//
Array.prototype.dim = function() {
	this.nrows = this.length;
	this.ncols = this[0].length;
};
//
// Print an array in matrix format using the trace operator
//
Array.prototype.print = function() {
	if (this.nrows>1) {
		for (var i = 0; i < this.nrows; i++) {
			trace("["+this[i]+"]");
		}
	} else {
		trace("["+this+"]");
	}
	trace(" ");
};
//
// Transpose a matrix 
//
Array.prototype.transpose = function() {
	var output = Matrix(this.ncols, this.nrows);
	output.dim();
	if (this.nrows>1) {
		for (var i = 0; i < this.nrows; i++) {
			for (var j = 0; j < this.ncols; j++) {
				output[j][i] = this[i][j];
			}
		}
	} else {
		for (var i = 0; i < this.nrows; i++) {
			for (var j = 0; j < this.ncols; j++) {
				output[j] = this[i][j];
			}
		}
	}
	output.nrows = this.ncols;
	output.ncols = this.nrows;
	return (output);
};
//
// Create a square identity matrix based on dimensions of this matrix
//
Array.prototype.identity = function() {
	var output = Matrix(this.ncols, this.nrows);
	if (this.nrows == this.ncols) {
		for (var i = 0; i < this.nrows; i++) {
			for (var j = 0; j < this.ncols; j++) {
				if (i == j) {
					output[i][j] = 1;
				} else {
					output[i][j] = 0;
				}
			}
		}
		output.dim();
		return (output);
	}
};
//
//
// -------------------------------------------
// Add, Subtract, and Multiply a scalar by a matrix
// -------------------------------------------
//
// Add a scalar to a matrix
//
Array.prototype.scalarAdd = function(k) {
	var output = Matrix(this.ncols, this.nrows);
	for (var i = 0; i  <  this.nrows; i++) {
		for (var j = 0; j  <  this.ncols; j++) {
			output[i][j] = this[i][j]+k;
		}
	}
	output.dim();
	return (output);
};
//
// Subtract a scalar from a matrix
//
Array.prototype.scalarSubtract = function(k) {
	var output = Matrix(this.ncols, this.nrows);
	for (var i = 0; i  <  this.nrows; i++) {
		for (var j = 0; j  <  this.ncols; j++) {
			output[i][j] = this[i][j]-k;
		}
	}
	output.dim();
	return (output);
};
//
// Multiply a matrix by a scalar, k
//
Array.prototype.scalarMult = function(k) {
	var output = Matrix(this.ncols, this.nrows);
	for (var i = 0; i  <  this.nrows; i++) {
		for (var j = 0; j  <  this.ncols; j++) {
			output[i][j] = this[i][j]*k;
		}
	}
	output.dim();
	return (output);
};
//

[chris23]

I need to place this in separate posts. Sorry. continues below.
--Chris

// -------------------------------------------
// Add, Subtract, and Multiply two matrices
// -------------------------------------------
//
// Add two matrices together
//
Array.prototype.matrixAdd = function(matrixb) {
	if ((this.nrows == matrixb.nrows) && (this.ncols == matrixb.ncols)) {
		var output = Matrix(this.nrows, this.ncols);
		for (var i = 0; i  <  this.nrows; i++) {
			for (var j = 0; j  <  this.ncols; j++) {
				output[i][j] = this[i][j]+matrixb[i][j];
			}
		}
		output.dim();
		return (output);
	}
};
//
// Subtract two matrices 
//   this.matrix - matrixb
//
Array.prototype.matrixSubtract = function(matrixb) {
	if ((this.nrows == matrixb.nrows) && (this.ncols == matrixb.ncols)) {
		var output = Matrix(this.nrows, this.ncols);
		for (var i = 0; i  <  this.nrows; i++) {
			for (var j = 0; j  <  this.ncols; j++) {
				output[i][j] = this[i][j]-matrixb[i][j];
			}
		}
		output.dim();
		return (output);
	}
};
//
// Matrix Multiplication
// Multiply two matrices | this matrix is *post multiplied* by matrixb
// so, A.matrixMult(B) is equivalent to AB
//
Array.prototype.matrixMult = function(matrixb) {
	if (this.ncols == matrixb.nrows) {
		var output = Matrix(this.nrows, matrixb.ncols);
		for (var i = 0; i  <  this.nrows; i++) {
			for (var k = 0; k  <  matrixb.ncols; k++) {
				var sumproducts = 0;
				for (var j = 0; j  <  matrixb.nrows; j++) {
					if (this.ncols == 1) {
						sumproducts = sumproducts+(this[i]*matrixb[j][k]);
					}
					if (matrixb.ncols == 1) {
						sumproducts = sumproducts+(this[i][j]*matrixb[j]);
					}
					if ((this.ncols != 1) && (matrixb.ncols != 1)) {
						sumproducts = sumproducts+(this[i][j]*matrixb[j][k]);
					}
				}
				if (this.nrows == 1) {
					output[i] = sumproducts;
				} else {
					output[i][k] = sumproducts;
				}
			}
		}
		output.dim();
		return (output);
	}
};
//
// Matrix * Vector Multiplication
// Post-multiply a matrix by a vector 
// so, A.matrixVectorMult(v) is equivalent to Av
//
Array.prototype.matrixVectorMult = function(vector) {
	var kk = vector.length;
	var output = new Array();
		for (var k = 0; k  <  this.ncols; k++) {
			var sumproducts = 0;
			for (var j = 0; j  <  kk; j++) {
				sumproducts = sumproducts+(this[k][j]*vector[j]);
			}
			output[k] = sumproducts;
		}
	output.dim();
	return (output);
};

[chris23]

Here is the code for the matrix inverse. --Chris

// Estimate the inverse and determinant of a square matrix
//
// based on the pivot method as explained by 
// Carroll, Green, & Chaturvedi (1997). Mathematical tools for applied multivariate
// analysis.  San Diego: Academic Press. (pp. 66 - 68; 184-187)
//
// I didn't spend the time to make this code efficient or polished.  
// It does work, however!  The matrix in question cannot have a 
// 0 in the leading entry.
//
Array.prototype.Inverse = function() {
	// find dimensions of matrix
	A = this;
	A.nrows = A.length;
	A.ncols = A[0].length;
	var k = A.nrows;
	pivots = new Array(k);
	//
	// Concatenate this matrix with an identity matrix of same dimensions.
	// This new matrix will be called AI.
	var B = A.identity();
	AI = Matrix(k, (k*2));
	for (var i = 0; i  <  k; i++) {
		temp = A[i].concat(B[i]);
		for (var j = 0; j  <  (k*2); j++) {
			AI[i][j] = temp[j];
		}
	}
	// Copy AI k times.  The different copies will be used to save
	// the operations at different points in the pivot process
	// AI0 will contain the original version of AI
	AInull = Matrix(k, (k*2));
	for (var i = 1; i  <  (k+1); i++) {
		_root["AI"+i] = Matrix(k, (k*2));
	}
	AI0 = AI;
	//
	// Loop through each pivoting level
	//
	for (var i = 0; i  <  (k-1); i++) {
		// select pivot
		pivot = _root["AI"+i][i][i];
		pivots[i] = pivot;
		// get target row from original matrix
		targetRow = new Array();
		for (var jj = 0; jj  <  (k*2); jj++) {
			targetRow[jj] = _root["AI"+i][i][jj];
		}
		// divide this row by pivot to create targetRow
		for (var jj = 0; jj  <  (k*2); jj++) {
			targetRow[jj] = targetRow[jj]/pivot;
		}
		// assign target row to next matrix
		_root["AI"+(i+1)][i] = targetRow;
		//
		// Loop through remaing rows of this matrix level and
		// perform row operations
		for (var j = (i+1); j  <  k; j++) {
			temp1 = new Array(k*2);
			// select multiplier (first non-zero element of next row of previous matrix)
			multiplier = _root["AI"+i][j][i];
			// multiply multiplier by targetRow 
			for (var jj = 0; jj  <  (k*2); jj++) {
				temp1[jj] = targetRow[jj]*multiplier;
			}
			// fill in counterpart 
			var trow2 = new Array();
			for (var jj = 0; jj  <  (k*2); jj++) {
				trow2[jj] = _root["AI"+i][j][jj];
			}
			var temp2 = new Array(k*2);
			// subtract temp1 from counterpart row 
			for (var jj = 0; jj  <  (k*2); jj++) {
				temp2[jj] = trow2[jj]-temp1[jj];
			}
			// fill in counterpart
			_root["AI"+(i+1)][j] = temp2;
		}
		//end j
	}
	// end i
	//
	// complete the last pivot
	pivots[(k-1)] = _root["AI"+(k-1)][(k-1)][(k-1)];
	// Compute the determinant of the matrix as the product of the pivots
	det = 1;
	for (var jj = 0; jj  <  k; jj++) {
		det = det*pivots[jj];
	}
	//
	// Use backsubstitution to solve for the inverse
	// result will be saved in matrix "inverse"
	//
	inverse = Matrix(k, k);
	inverseCol = 0;
	// find first row of backwards inverse matrix
	targetRow = new Array();
	for (var jj = 0; jj  <  (k*2); jj++) {
		targetRow[jj] = _root["AI"+(k-1)][k-1][jj];
	}
	// divide targetRow by pivots[k]
	for (var jj = 0; jj  <  (k*2); jj++) {
		targetRow[jj] = targetRow[jj]/pivots[k-1];
	}
	// insert targetRow into AI level k+1
	for (var jj = 0; jj  <  (k*2); jj++) {
		_root["AI"+(k)][0][jj] = targetRow[jj];
	}
	// loop through each column of the final matrix
	//
	for (var j = k; j  <  (k+k-0); j++) {
		endvec = new Array(k-1);
		// fill endvec with zeros
		for (var jj = 0; jj  <  (k-1); jj++) {
			endvec[jj] = 0;
		}
		// allow the last element of endvec to equal last piece
		endvec[k-1] = _root["AI"+(k)][0][j];
		//
		// loop thru and calculate inverse within each column
		for (var i = (k-2); i>=0; i--) {
			// pull out row i from AI level i+1
			temp = new Array();
			temp = _root["AI"+(i+1)][i];
			// multiply each element of that row by each element of endvec. sum products
			tempsum = 0;
			for (var jj = 0; jj  <  k; jj++) {
				tempsum = tempsum+_root["AI"+(i+1)][i][jj]*endvec[jj];
			}
			// subject tempsum from temp
			endvec[i] = _root["AI"+(i+1)][i][j]-tempsum;
		}
		// save result to correct part of inverse matrix
		for (var jj = 0; jj  <  k; jj++) {
			inverse[jj][inverseCol] = endvec[jj];
		}
		inverseCol = inverseCol+1;
	}
	// end j
	trace("determinant = "+det);
	return (inverse);
};

[eyezberg]

Welcome on bo(a)rd, Einstein!
I'm impressed, do you like Math?

[chris23]

Yeah, math is boring, isn't it?

I'm not a big fan of mathematics for mathematics sake, but I use a lot of math in my work. Sigh.

[Jubba]

what exactly is this useful for? i am not familiar with matrices...

[chris23]

Well, matrix operations are useful for a number of applications. As a researcher, the most important application or me is the analysis of multi-variable data sets. A substantial portion of contemporary statistical techniques (such as multiple regression, factor analysis, linear structural modeling) are built on the mathematics of matrices. In order to use Flash as an interactive medium for on-line statistical applications, it is necessary to be able to manipulate matrices easily and to calculate the inverse of a matrix.

For someone interested in graphics and interactive design (i.e., the non-nerd ), matrices may be used to scale, rotate, and transform information (e.g., the x and y coordinates of movie clips or the scale of the clip). Many of the 3D transformations that are used in Flash and other applications are based on matrices (or on mathematics that could be manipulated more efficiently in matrix form). For example, a matrix containing the paired x and y coordinates of many movie clips could be passed through a single transformation matrix that would yield the rotated coordinates all of the clips (in one command) through space. I’ll try to create a concrete example later this week to illustrate.

[pom]

I haven't read all of it but it looks very interesting, chris

Concerning matrices, Brandon Williams did a very good job trying to explain how to use them on his site (link in Best of Kirupa) and in his tutes, for example this one: http://www.ultrashock.com/tutorials/..._tutorial.html

pom

[chris23]

That's a great overview of matrix algebra, Pom! Thanks for posting that link.

Later this week I'll write a mini-overview of how matrix algebra can be used for some simple statistical operations. Just for fun, however, I started applying some of the matrix functions to do 3D stuff. Since the other tutorial explains the matrix operations in depth, I won't elaborate too much here except to highlight some of the ways in which using matrices can be handy.

If we have a single movie clip and want to rotate it in 3D space, we can multiply the x, y, and z coordinates by the appropriate trig functions. For example, if we wanted to rotate a clip around the z axis, we would apply the following mathematics to the clip's coordinates:

x' = cos(r)*x - sin(r)*y
y' = sin(r)*x + cos(r)*y
z' = z

Although this is a simple example and doesn't necessarily benefit from the use of matrices, it is possible to represent this set of operations in matrix form:

V' = AV

where V' is a 3 x 1 matrix containing the transformed coordinates, A is a 3 x 3 matrix containing the transformation weights/functions, and V is the 3 x 1 matrix of original coordinates:

V'   =         A             V
[x'] = [cos(r)  -sin(r)  0] [x]
[y'] = [sin(r)   cos(r)  0] [x]
[z'] = [0         0      1] [x].

Once we specify the general matrix structure, we can apply the same equation to rotate multiple clips at once. All that needs to be done is expand matrix V to include a different column containing the coordinates of each clip.

     V'           =         A                 V
[x1' x2' x3' x4'] = [cos(r)  -sin(r)  0] [x1 x2 x3 x4]
[y1' y2' y3' y4'] = [sin(r)   cos(r)  0] [y1 y2 y3 y4]
[z1' z2' z3' z4'] = [0         0      1] [z1 z2 z3 z4].

In other words, all of the coordinate information for the various movie clips can be placed in a single matrix, and all of the coordinates (or just selected coordinates, if desired) can be transformed in a simple operation.

Another advantage of using matrices in this context is that it provides an easy way to combine rotations. In other words, if one wants to rotate a set of clips along the x and y axis simultaneously, one can easily add (or subtract, depending on the direction of rotation) the two transformation matrices together using some of these new functions. If one wanted to stretch of shear the space, one could multiply the different transformation matrices--allowing for really kewl forms of motion.

As I noted above, I won't beat this topic to death since there are already good tutorials and sources on using 3D in ActionScript. I put together a simple swf file, however, that illustrates some of the rotations. This file (see link below) uses the matrix functions described previously, and all clips are rotated simultaneously. To rotate the points, simply click on the appropriate button. Click repeatedly to advance the rotation more quickly. The last rotation is based on multiplying--using matrix multiplication--the transformation matrix for x and y rotation. The effect is pretty neat--it kinds of reminds me of leaves blowing about in an alley way.

Just for fun, I'll post some matrix stuff on statistics in a few days.
Chris

link to file:
http://tigger.uic.edu/~fraley/matrix-rotate.htm