# The Matrix Isn't Real

Last time we saw that because complex numbers support all of the familiar arithmetic operations it was possible to define arithmetic operations for vectors having complex rather than real elements[1], although for the sake of one's sanity it's probably best not to think too much about the directions in which they point.
Given that our implementation of them with the ak.complexVector class simply required us to use our own overloaded arithmetic operators rather than JavaScript's native ones, it's not unreasonable to expect that doing the same thing for matrices should be no more difficult.

### A Brief Refresher

When we first covered matrices[2] I described them as representing linear transformations of vectors that result in vectors whose elements are simply linear combinations of those of the original vector. For example
$\begin{pmatrix} a & b\\ c & d \end{pmatrix} \times \begin{pmatrix} x\\ y \end{pmatrix} = \begin{pmatrix} ax + by\\ cx + dy \end{pmatrix}$
and
$\begin{pmatrix} x & y \end{pmatrix} \times \begin{pmatrix} a & b\\ c & d \end{pmatrix} = \begin{pmatrix} ax + cy\\ bx + dy \end{pmatrix}$
Note that writing a vector that appears on the right hand side of a multiplication with a matrix in a column and one that appears on the left hand side in a row is purely a notational convention and the vectors
$\begin{pmatrix} x\\ y \end{pmatrix}$
and
$\begin{pmatrix} x & y \end{pmatrix}$
are identical. That such a convention is used is because pre-multiplying $$n$$ dimensional vectors and $$n$$ row by one column matrices by $$m$$ row by $$n$$ column matrices are extremely similar operations, as are post-multiplying $$n$$ dimensional vectors and one row by $$n$$ column matrices by $$n$$ row by $$m$$ column matrices.

If we label the element in the $$i^{\mathrm{th}}$$ row and the $$j^{\mathrm{th}}$$ column of a matrix $$\mathbf{M}$$ with $$m_{ij}$$ and the $$i^{\mathrm{th}}$$ element of a vector $$\mathbf{v}$$ with $$v_i$$, then the rules for multiplying matrices and vectors are given by
\begin{align*} (\mathbf{M} \times \mathbf{v})_i &= \sum_j m_{ij} \times v_j\\ (\mathbf{v} \times \mathbf{M})_j &= \sum_i v_i \times m_{ij}\\ (\mathbf{M} \times \mathbf{N})_{ij} &= \sum_k m_{ik} \times n_{kj} \end{align*}
where $$\sum$$ is the summation sign.
Note that this implies that if a vector is pre-multiplied by a matrix then it must have the same number of dimensions as the matrix has columns and that if one is post-multiplied by a matrix then the number of its dimensions must equal the number of the matrix's rows. Furthermore, if two matrices are multiplied then the one on the left hand side must have the same number of columns as the one on the right hand side has rows.

If compare post-multiplication of a vector by a matrix with that of a matrix with one row, let's call it $$\mathbf{V}$$, we can see why the notational convention is an appealing one
\begin{align*} (\mathbf{v} \times \mathbf{M})_j &= \sum_i v_i \times m_{ij}\\ (\mathbf{V} \times \mathbf{M})_{1j} &= \sum_i v_{1i} \times m_{ij} \end{align*}
Now as a refresher in matrix arithmetic this is admittedly rather inadequate, but hopefully this purely algebraic definition of multiplication has put out of your mind any philosophical misgivings that you might have had about the arithmetic of matrices with complex elements.

### A Complex Matrix Class

Just as our ak.complexVector class had methods in common with both the ak.vector and ak.complex classes, so our complex matrix class, ak.complexMatrix, has methods in common with both the ak.matrix and ak.complex classes, as illustrated by listing 1.

Listing 1: ak.complexMatrix
ak.COMPLEX_MATRIX_T = 'ak.complexMatrix';

function ComplexMatrix(){}
ComplexMatrix.prototype = {
TYPE: ak.COMPLEX_MATRIX_T, valueOf: function(){return ak.NaN;}
};

ak.complexMatrix = function() {
var m     = new ComplexMatrix();
var state = {rows: 0, cols: 0, elems: []};
var arg0  = arguments[0];

constructors[ak.nativeType(arg0)](state, arg0, arguments);

m.rows = function() {return state.rows;};
m.cols = function() {return state.cols;};

m.at = function(r, c) {return state.elems[Number(r)*state.cols+Number(c)];};

m.re = function() {
return ak.matrix(state.rows,
state.cols,
state.elems.map(function(x){return x.re();}));
};

m.im = function() {
return ak.matrix(state.rows,
state.cols,
state.elems.map(function(x){return x.im();}));
};

m.toExponential = function(d) {
state.cols,
state.elems,
function(x){return x.toExponential(d);});
};

m.toFixed = function(d) {
state.cols,
state.elems,
function(x){return x.toFixed(d);});
};

m.toPrecision = function(d) {
state.cols,
state.elems,
function(x){return x.toPrecision(d);});
};

return Object.freeze(m);
};

var constructors = {};


As usual, we have an internal state object that is initialised by dispatching to a constructors object, in this case to set the rows, cols and elems, with the latter being a one dimensional array of the matrix elements arranged in a row-wise fashion.
Furthermore, we have the property access methods rows, cols and at to provide access to the number of rows, the number of columns and the elements of the complex matrix, and the methods re and im to extract matrices containing the real and imaginary parts of its elements respectively. Note how the at method converts the row and column into an offset into the elems.
Finally, we have the type conversion methods toArray, toString, toExponential, toFixed and toPrecision to convert the matrix into a native JavaScript array of arrays and string representations with the standard number formats and which depend upon the helper functions toArray and toString given in listing 2.

Listing 2: Type Conversion Helper Functions
function toArray(rows, cols, elems) {
var a = new Array(rows);
var off = 0;
var r, row, c;

for(r=0;r<rows;++r) {
row = new Array(cols);
for(c=0;c<cols;++c) row[c] = elems[off++];
a[r] = row;
}
return a;
}

function toString(rows, cols, elems, f) {
var s = [];
var off = 0;
var r, c;

s.push('[');

for(r=0;r<rows;++r) {
s.push('[');

for(c=0;c<cols;++c) {
s.push(f ? f(elems[off++]) : elems[off++]);
s.push(',');
}

s.pop();
s.push(']');
s.push(',');
}

s.pop();
s.push(']');

return s.join('');
}


Here we convert the one dimensional elems into two dimensional array and string representations by keeping track of the offset into it with off as we iterate over the rows and columns of the matrix.

### ak.complexMatrix Constructors

We shall again follow the approach that we did for the ak.complexVector class and support construction of ak.complexMatrix objects in the same ways that we do for the ak.matrix class, albeit allowing the separate specification of the real and imaginary parts of its elements as well as that of their complex values.

The simplest way to initialise an ak.complexMatrix object is with an array of arrays of the ak.complex values in its rows or, failing that, with a pair of arrays of arrays containing the real and imaginary parts of the elements in its rows respectively. As we did with ak.matrix, we assume that the number of columns is equal to the number of elements of the first array contained in those arrays, as illustrated by listing 3.

Listing 3: ak.complexMatrix Array Constructors
constructors[ak.ARRAY_T] = function(state, arr, args) {
var arg1 = args[1];

state.rows = arr.length;
state.cols = state.rows ? arr[0].length : 0;
state.elems = new Array(state.rows*state.cols);

constructors[ak.ARRAY_T][ak.nativeType(arg1)](state, arr, arg1);
};

constructors[ak.ARRAY_T][ak.UNDEFINED_T] = function(state, arr) {
var off = 0;
var r, row, c;

for(r=0;r<state.rows;++r) {
row = arr[r];
for(c=0;c<state.cols;++c) state.elems[off++] = ak.complex(row[c]);
}
};

constructors[ak.ARRAY_T][ak.ARRAY_T] = function(state, re, im) {
var off = 0;
var r, rer, imr, c;

if(im.length!==state.rows || state.rows>0 && im[0].length!==state.cols) {
throw new Error('dimensions mismatch in ak.complexMatrix');
}

for(r=0;r<state.rows;++r) {
rer = re[r];
imr = im[r];
for(c=0;c<state.cols;++c) {
state.elems[off++] = ak.complex(Number(rer[c]), Number(imr[c]));
}
}
};


As was the case with ak.matrix, we perform no checks that the arrays are all the same length, ignoring elements in row arrays that are longer than the first and reading undefined values from past the end of those that are shorter. However, since ak.complex objects cannot be initialised with undefined values the array constructors will throw an exception if the argument contains arrays that are shorter than the first rather than populate them with NaNs, as did the ak.matrix array constructor and as will the real and imaginary arrays constructor.

Like we did for ak.complexVector, we shall take the same approach when initialising with the number of rows and columns and a value or values to set the elements to, as illustrated by listing 4.

Listing 4: ak.complexMatrix Size Constructors
constructors[ak.NUMBER_T] = function(state, rows, args) {
var arg1 = args[1];

state.rows = rows;
constructors[ak.NUMBER_T][ak.nativeType(arg1)](state, arg1, args);
};

constructors[ak.NUMBER_T][ak.NUMBER_T] = function(state, cols, args) {
var arg2 = args[2];

state.cols = cols;
state.elems = new Array(state.rows*state.cols);
constructors[ak.NUMBER_T][ak.NUMBER_T][ak.nativeType(arg2)](state, arg2, args);
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.ARRAY_T] = function(state, vals, args) {
var arg3 = args[3];
constructors[ak.NUMBER_T][ak.NUMBER_T][ak.ARRAY_T][ak.nativeType(arg3)](state, vals, arg3);
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.ARRAY_T][ak.ARRAY_T] = function(state, re, im) {
var n = state.rows*state.cols;
var i;

for(i=0;i<n;++i) state.elems[i] = ak.complex(Number(re[i]), Number(im[i]));
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.ARRAY_T][ak.UNDEFINED_T] = function(state, vals) {
var n = state.rows*state.cols;
var i;

for(i=0;i<n;++i) state.elems[i] = ak.complex(vals[i]);
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.FUNCTION_T] = function(state, f) {
var off = 0;
var r, c;

for(r=0;r<state.rows;++r) {
for(c=0;c<state.cols;++c) {
state.elems[off++] = ak.complex(f(r, c));
}
}
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T] = function(state, val, args) {
var arg3 = args[3];
constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.nativeType(arg3)](state, val, arg3);
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T] = function(state, re, im) {
var n = state.rows*state.cols;
var c = ak.complex(re, im);
var i;

for(i=0;i<n;++i) state.elems[i] = c;
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.NUMBER_T][ak.UNDEFINED_T] = function(state, val) {
var n = state.rows*state.cols;
var i;

val = ak.complex(val);
for(i=0;i<n;++i) state.elems[i] = val;
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.OBJECT_T] = function(state, obj) {
var n = state.rows*state.cols;
var i;

obj = ak.complex(obj);
for(i=0;i<n;++i) state.elems[i] = obj;
};

constructors[ak.NUMBER_T][ak.NUMBER_T][ak.UNDEFINED_T] = function(state) {
var n = state.rows*state.cols;
var c = ak.complex(0);
var i;

for(i=0;i<n;++i) state.elems[i] = c;
};


Note that included in these is the positional constructor in which the elements are set to the results of a function of their row and column indices. Furthermore, if the number of rows and columns are followed by an array, or arrays, of values they are expected to have the same row-wise ordering that the elems member of the internal state object has.

It is also convenient to initialise ak.complexMatrix object with other ak.complexMatrix objects, or at least with objects having similar interfaces, or with pairs of ak.matrix objects, or at least with pairs of objects having similar interfaces, representing the real and imaginary parts, as enabled by the constructors given by listing 5.

Listing 5: ak.complexMatrix Object Constructors
constructors[ak.OBJECT_T] = function(state, obj, args) {
var arg1 = args[1];
constructors[ak.OBJECT_T][ak.nativeType(arg1)](state, obj, arg1);
}

constructors[ak.OBJECT_T][ak.OBJECT_T] = function(state, re, im) {
var rrows = re.rows;
var rcols = re.cols;
var irows = im.rows;
var icols = im.cols;
var off = 0;
var elems, r, c;

rrows = (ak.nativeType(rrows)===ak.FUNCTION_T) ? Number(rrows()) : Number(rrows);
rcols = (ak.nativeType(rcols)===ak.FUNCTION_T) ? Number(rcols()) : Number(rcols);

irows = (ak.nativeType(irows)===ak.FUNCTION_T) ? Number(irows()) : Number(irows);
icols = (ak.nativeType(icols)===ak.FUNCTION_T) ? Number(icols()) : Number(icols);

if(irows!==rrows || icols!==rcols) {
throw new Error('dimensions mismatch in ak.complexMatrix');
}

elems = new Array(rrows*rcols);

for(r=0;r<rrows;++r) {
for(c=0;c<rcols;++c) {
elems[off++] = ak.complex(Number(re.at(r, c)), Number(im.at(r, c)));
}
}

state.rows  = rrows;
state.cols  = rcols;
state.elems = elems;
};

constructors[ak.OBJECT_T][ak.UNDEFINED_T] = function(state, obj) {
var rows = obj.rows;
var cols = obj.cols;
var off = 0;
var elems, r, c;

rows = (ak.nativeType(rows)===ak.FUNCTION_T) ? Number(rows()) : Number(rows);
cols = (ak.nativeType(cols)===ak.FUNCTION_T) ? Number(cols()) : Number(cols);

elems = new Array(rows*cols);

for(r=0;r<rows;++r) {
for(c=0;c<cols;++c) {
elems[off++] = ak.complex(obj.at(r, c));
}
}

state.rows  = rows;
state.cols  = cols;
state.elems = elems;
};


Finally, we shall allow the explicit construction of complex identity and diagonal matrices in the same way that we did for normal matrices; by passing a string identifying the type of matrix as the first argument, as shown in listing 6.

Listing 6: ak.complexMatrix Special Form Constructors
constructors[ak.STRING_T] = function(state, type, args) {
constructors[ak.STRING_T][type.toUpperCase()](state, args);
};

constructors[ak.STRING_T]['IDENTITY'] = function(state, args) {
var arg1 = args[1];
constructors[ak.STRING_T]['IDENTITY'][ak.nativeType(arg1)](state, arg1);
};

constructors[ak.STRING_T]['IDENTITY'][ak.NUMBER_T] = function(state, n) {
var elems = new Array(n*n);
var off = 0;
var zero = ak.complex(0);
var one = ak.complex(1);
var r, c;

for(r=0;r<n;++r) {
for(c=0;c<r;++c)   elems[off++] = zero;
elems[off++] = one;
for(c=r+1;c<n;++c) elems[off++] = zero;
}

state.rows  = n;
state.cols  = n;
state.elems = elems;
};

constructors[ak.STRING_T]['DIAGONAL'] = function(state, args) {
var a = new Array(args.length-1);
var off = 0;
var zero = ak.complex(0);
var i, v, r, c, n, elems;

for(i=1;i<args.length;++i) a[i-1] = args[i];
v = ak.complexVector.apply(null, a);
n = v.dims();
elems = new Array(n*n);

for(r=0;r<n;++r) {
for(c=0;c<r;++c)   elems[off++] = zero;
elems[off++] = v.at(r);
for(c=r+1;c<n;++c) elems[off++] = zero;
}

state.rows  = n;
state.cols  = n;
state.elems = elems;
};


Identity matrices, being those that have elements equal to one when their row indices equals their column indices and equal to zero otherwise, must be square so we only need a single number to specify the number of both the rows and the columns.
For diagonal matrices, being those whose elements with different row and column indices are equal to zero and which also must be square, we populate an array with all of the arguments except the first, being the one that indicated that we wanted to construct a diagonal matrix in the first place, and use it as the arguments for the construction of an ak.complexVector object which we then use to initialise the diagonal elements. This allows us to use a similar form of construction for diagonal complex matrices that we have for complex vectors.

To implement complex matrix arithmetic we need simply replace JavaScript's native arithmetic operators with our own overloaded ones in the code from ak.matrix, much as we did for ak.complexVector. By way of an example, listing 7 show how we can do this for multiplication, the rules of which we covered above.

Listing 7: ak.complexMatrix Multiplication
function mul(w0, w1) {
var r0 = w0.rows();
var c0 = w0.cols();
var r1 = w1.rows();
var c1 = w1.cols();
var a = new Array(r0*c1);
var off = 0;
var i, j, k, s;

if(c0!==r1) throw new Error('dimensions mismatch in ak.complexMatrix mul');

for(i=0;i<r0;++i) {
for(j=0;j<c1;++j) {
s = ak.complex(0);
for(k=0;k<c0;++k) s = ak.add(s, ak.mul(w0.at(i, k), w1.at(k, j)));
a[off++] = s;
}
}
return ak.complexMatrix(r0, c1, a);
}

function mulWZ(w, z) {
var rows = w.rows();
var cols = w.cols();
var dims = z.dims();
var a = new Array(rows);
var r, c, s;

if(dims!==cols) throw new Error('dimensions mismatch in ak.complexMatrix mul');

for(r=0;r<rows;++r) {
s = 0;
for(c=0;c<cols;++c) s = ak.add(s, ak.mul(w.at(r, c), z.at(c)));
a[r] = s;
}
return ak.complexVector(a);
}

function mulZW(z, w) {
var rows = w.rows();
var cols = w.cols();
var dims = z.dims();
var a = new Array(cols);
var r, c, s;

if(dims!==rows) throw new Error('dimensions mismatch in ak.complexMatrix mul');

for(c=0;c<cols;++c) {
s = 0;
for(r=0;r<rows;++r) s = ak.add(s, ak.mul(z.at(r), w.at(r, c)));
a[c] = s;
}
return ak.complexVector(a);
}



Note that since our overloaded operators work for mixed real and complex arithmetic, we can reuse the implementations of multiplication for the various combinations of real and complex matrices and vectors, as demonstrated in part by program 1.

Program 1: Complex Matrix Multiplication

For a more complicated example, we can implement the inverse of a complex matrix in exactly the same way, as shown by listing 8.

Listing 8: ak.complexMatrix Inverse
function inv(w) {
var n = w.rows();
var a = w.toArray();
var i = ak.complexMatrix('identity', n).toArray();
var r1, r2, c, arow1, irow1, abs1, arow2, irow2, abs2, z;

if(w.cols()!==n) {
throw new Error('non-square matrix in ak.complexMatrix inv');
}

for(r1=0;r1<n;++r1) {
arow1 = a[r1];
irow1 = i[r1];
abs1  = ak.abs(arow1[r1]);

for(r2=r1+1;r2<n;++r2) {
arow2 = a[r2];
abs2  = ak.abs(arow2[r1]);

if(abs2>abs1) {
irow2 = i[r2];

a[r2] = arow1; arow1 = arow2; a[r1] = arow1;
i[r2] = irow1; irow1 = irow2; i[r1] = irow1;
abs1  = abs2;
}
}

if(abs1===0) {
throw new Error('singular matrix in ak.complexMatrix inv');
}
z = arow1[r1];

for(c=r1+1;c<n;++c) arow1[c] = ak.div(arow1[c], z);
for(c=0;   c<n;++c) irow1[c] = ak.div(irow1[c], z);

for(r2=0;r2<n;++r2) {
if(r2!==r1) {
arow2 = a[r2];
irow2 = i[r2];
z = arow2[r1];

for(c=r1+1;c<n;++c) arow2[c] = ak.sub(arow2[c], ak.mul(z, arow1[c]));
for(c=0;   c<n;++c) irow2[c] = ak.sub(irow2[c], ak.mul(z, irow1[c]));
}
}
}
return ak.complexMatrix(i);
}



Whilst you will no doubt recall that the inverse of a square matrix $$\mathbf{M}$$, denoted by $$\mathbf{M}^{-1}$$, satisfies
$\mathbf{M}^{-1} \times \mathbf{M} = \mathbf{M} \times \mathbf{M}^{-1} = \mathbf{I}$
where $$\mathbf{I}$$ is the identity matrix, you may not remember the details of the Gaussian elimination algorithm for finding it, of which this is an example. Since I have absolutely no intention of going over it again, if you want a description of how it works I suggest you refer back to our treatment of real matrices[2]. However, that it works even for complex matrices is demonstrated by program 2.

Program 2: The Complex Matrix Inverse

As was the case for complex vectors, we need to be a little careful when it comes to the magnitude of complex matrices. Specifically, rather than take the square root of the sum of the squares of the elements, which for complex values can be zero for non-zero values, we must take the square root of the sum of the squares of their magnitudes, which can't
$\sqrt{\sum_i \sum_j \left|m_{ij}\right|^2}$
where $$m_{ij}$$ is the element in the $$i^\mathrm{th}$$ row and the $$j^\mathrm{th}$$ column, as is done in the overload of ak.abs shown in listing 9.

Listing 9: ak.complexMatrix Magnitude
function abs(w) {
var s = 0;
var rows = w.rows();
var cols = w.cols();
var r, c, x;

for(r=0;r<rows;++r) {
for(c=0;c<cols;++c) {
s += Math.pow(ak.abs(w.at(r, c)), 2);
}
}
return Math.sqrt(s);
}



Since these matrices are complex valued they also have conjugates, in which the imaginary parts are negated, as implemented in listing 10.

Listing 10: ak.complexMatrix And ak.matrix Conjugates
function conj(w) {
var rows = w.rows();
var cols = w.cols();
var a = new Array(rows*cols);
var off = 0;
var r, c;

for(r=0;r<rows;++r) {
for(c=0;c<cols;++c) {
a[off++] = ak.conj(w.at(r, c));
}
}
return ak.complexMatrix(rows, cols, a);
}

function conjM(m) {
return m;
}



Note that we also implement the conjugate for real matrices, which is simply the matrix itself since it has no imaginary part.

In addition to the transpose, in which a matrix is flipped over to swap its rows with its columns, a complex matrix has a Hermitian conjugate or adjoint, often written as $$\mathbf{M}^\mathrm{H}$$, which is the conjugate of its transpose. Some important properties of symmetric real matrices, which are equal to their transposes, are also to be found in Hermitian complex matrices, which are equal to their adjoints.
Whilst I shan't yet specify what those properties actually are, I shall assert that they are sufficiently important to justify adding an adjoint function to our overloaded arithmetic operators, as is done in listing 11.

Listing 11: ak.complexMatrix And ak.matrix Adjoints
if(!ak.adjoint) ak.adjoint = function(x) {return ak.adjoint[ak.type(x)](x)};

var rows = w.rows();
var cols = w.cols();
var a = new Array(rows*cols);
var off = 0;
var r, c;

for(r=0;r<rows;++r) {
for(c=0;c<cols;++c) {
a[off++] = ak.conj(w.at(c, r));
}
}
return ak.complexMatrix(cols, rows, a);
}



Having implemented the adjoint for complex matrices, we also do so for real matrices too for which it is trivially the same as the transpose.

The rest of the arithmetic operations for complex matrices simply replace JavaScript's arithmetic operators with our own, just like the overloads of ak.mul and ak.inv did, and so aren't really worth examining in detail. They are shown in listing 12, however.

var JACOBI_DECOMPOSITION_T = 'ak.jacobiDecomposition';



Note that, just as we did last time, we are using a local type definition for the Jacobi decomposition overloads[3] so that the user may choose whether or not to load the file in which ak.jacobiDecomposition is defined.

The implementations of all of the complex matrix overloads can be found in ComplexMatrix.js and you check out some of those that we haven't covered in program 3.

Program 3: Miscellaneous Complex Matrix Operators

Of course, the interesting aspects of linear algebra involve trying to find special forms for matrices that allow us to perform more sophisticated operations, much as we did with the Jacobi decomposition. In the next few posts we shall take a look at an admittedly much simpler one to see how easily we can extend them to complex matrices too.

### References

[1] What Are You Pointing At?, www.thusspakeak.com, 2015.

[2] The Matrix Recoded, www.thusspakeak.com, 2014.

[3] Conquering The Eigen, www.thusspakeak.com, 2014.

Hi Richard

Please could you elucidate the self adjoint operator of the Green function.

All the best Rupert