## Saturday, July 19, 2014

### C# Matrix Inversion With Latex Output

Matrix inversion is important in the domains of solving a large number of problems. Many types of problems can be expressed and then solved using matrices. In my areas of work, I most frequently encounter matrix inversion problems when working with regression problems and math problems involving affine transformations (ex. those used in cryptography applications).

Anyway, from my first college level linear algebra class, I typed up most of my math homework using Latex, and I worked out a reasonable way to show Matrix inversion step by step using row multiplications and row additions. Even though it's been a number of years since I've taken a course in linear algebra, I've wanted to build an application that outlines all of the single row multiplication and multiple row additions that are required to find the inverse for a matrix of arbitrary rank (rank is the number of diagonal elements in a square matrix).

Researchers might use this application to help them find a solution to a specific problem. Students might use the application to check their homework or to learn more about Gauss Jordan elimination.
Using this application, in a few posts I plan to answer the following questions:
What is the inverse of a 1x1 Matrix?
What is the inverse of a 2x2 Matrix?
What is the inverse of a 3x3 Matrix?
What is the inverse of a 4x4 Matrix?
Depending on the input matrix, initial row substitutions might be necessary to ensure that all diagonal elements are non-zero during the elimination process, since the application has no knowledge of the input matrix (this is directly entered into the generated latex code). Additionally, some steps may need to be removed from the output as they might be unnecessary (ex. initial zero elements off of the diagonal of the input matrix).

This application can, with sufficient memory, disk space, and compiled as 64 bit, invert large matrices, but Latex has an inherent limitation of paper sizes around 200 in wide. Testing would suggest that any matrix over rank 4 will probably overrun this limitation. Additionally, the disk size for all of the intermediate steps becomes really large (my machine tops out computing an 8x8, with an output file size > 10 GB). Latex also has some difficulty compiling these large files (numerous changes to texmf.cnf). These settings are about the highest that I can apply to my system:

main_memory    = 10000000
extra_mem_top  = 10000000
max_strings    = 2000000
hash_extra     = 10000000
pool_size      = 25000000

If you need to compile the .tex files linked to my blog or generated from my tex generating applications, you should probably use Tex Live: http://tug.org/texlive/

Here is an example of the tex output of the application for a 2x2 matrix (note that you can change the variables that appear in the final output):

\documentclass{article}

% This is the output from LatexMatrixInverse 1.0 for a matrix with rank 2
% http://mikemstech.blogspot.com

\usepackage{geometry}

% Note: you should probably use pdflatex to compiile this file.
% Other processors are known to have some issues with using
% 'geometry' to set paper size

% Adjust the page size here if output is wrapping in a bad way.
% Default is 8.5 x 11 in (Letter)
\geometry{papersize={16in,11in}}

%Import AMS Latex packages
\usepackage{amsmath, amssymb}
\setcounter{MaxMatrixCols}{5}

%Variable definition

\begin{document}
% Definition of initial A
% A row 1
\newcommand{\ARbCb}{a_{1,1}}
\newcommand{\ARbCc}{a_{1,2}}

% A row 2
\newcommand{\ARcCb}{a_{2,1}}
\newcommand{\ARcCc}{a_{2,2}}

% Definition of initial B
\newcommand{\BRb}{b_{1}}
\newcommand{\BRc}{b_{2}}

LatexMatrixInverse 1.0 Output for rank 2, ShowIntermediateSteps is True.

Given the following initial matrices:
\begin{equation*}
A = \begin{pmatrix}\ARbCb
&\ARbCc
\\\ARcCb
&\ARcCc
\end{pmatrix}B = \begin{pmatrix}\BRb\\ \BRc\end{pmatrix}\end{equation*}
We want to find $A^{-1}$ and $A^{-1} B$...
\begin{equation*}
\left ( \begin{array}{cc|cc|c}\ARbCb
&\ARbCc
&1
&0
&\BRb\\
\ARcCb
&\ARcCc
&0
&1
&\BRc\end{array} \right )
\end{equation*}
\begin{equation*}
\xrightarrow{\frac{1}{\ARbCb} R_{1}}
\left ( \begin{array}{cc|cc|c}1
&\frac{\ARbCc}{\ARbCb}
&\frac{1}{\ARbCb}
&0
&\frac{\BRb}{\ARbCb}\\
\ARcCb
&\ARcCc
&0
&1
&\BRc\end{array} \right )
\end{equation*}
\begin{equation*}
\xrightarrow{R_{2} - \left ( \ARcCb\right ) R_{1}}
\left ( \begin{array}{cc|cc|c}1
&\frac{\ARbCc}{\ARbCb}
&\frac{1}{\ARbCb}
&0
&\frac{\BRb}{\ARbCb}\\
0
&\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )
&0 - \left (\ARcCb \right ) \left ( \frac{1}{\ARbCb}\right )
&1
&\BRc - \left (\ARcCb \right ) \left ( \frac{\BRb}{\ARbCb}\right )\end{array} \right )
\end{equation*}
\begin{equation*}
\xrightarrow{\frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )} R_{2}}
\left ( \begin{array}{cc|cc|c}1
&\frac{\ARbCc}{\ARbCb}
&\frac{1}{\ARbCb}
&0
&\frac{\BRb}{\ARbCb}\\
0
&1
&\frac{0 - \left (\ARcCb \right ) \left ( \frac{1}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
&\frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
&\frac{\BRc - \left (\ARcCb \right ) \left ( \frac{\BRb}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\end{array} \right )
\end{equation*}
\begin{equation*}
\xrightarrow{R_{1} - \left ( \frac{\ARbCc}{\ARbCb}\right ) R_{2}}
\left ( \begin{array}{cc|cc|c}1
&0
&\frac{1}{\ARbCb} - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{0 - \left (\ARcCb \right )
\left ( \frac{1}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )
&0 - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )
&\frac{\BRb}{\ARbCb} - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{\BRc - \left (\ARcCb \right ) \left ( \frac{\BRb}{\ARbCb}\right )}{\ARcCc
- \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )\\
0
&1
&\frac{0 - \left (\ARcCb \right ) \left ( \frac{1}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
&\frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
&\frac{\BRc - \left (\ARcCb \right ) \left ( \frac{\BRb}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\end{array} \right )
\end{equation*}
\begin{equation*}
A^{-1} = \begin{pmatrix}\frac{1}{\ARbCb} - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{0 - \left (\ARcCb \right ) \left ( \frac{1}{\ARbCb}\right )}{\ARcCc
- \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )
&
0 - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )
\\\frac{0 - \left (\ARcCb \right ) \left ( \frac{1}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
&
\frac{1}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}
\end{pmatrix}A^{-1}B = \begin{pmatrix}\frac{\BRb}{\ARbCb} - \left (\frac{\ARbCc}{\ARbCb} \right ) \linebreak \left ( \frac{\BRc - \left (\ARcCb \right )
\left ( \frac{\BRb}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\right )\\ \frac{\BRc -
\left (\ARcCb \right ) \left ( \frac{\BRb}{\ARbCb}\right )}{\ARcCc - \left (\ARcCb \right ) \left ( \frac{\ARbCc}{\ARbCb}\right )}\end{pmatrix}
\end{equation*}
\end{document}


Below is a screenshot of the PDF output from the 2x2 file:

Obviously, the end user would need to perform any simplification at the end.

Below is the application in C#

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace LatexMatrixInverse
{
public class Program
{
public int rank = 0;

public string[,] A,
A_inverse;
public string[] B;

public string ArrayFormat = null;

public bool ShowIntermediateSteps = true;

public bool FinalOutputAsMatrix = true;

static void Main(string[] args)
{
Program p = new Program();

if (args.Length == 1)
{
p.rank = 0;

//arg 0 is rank
try { p.rank = int.Parse(args[0]); }
catch (Exception) { PrintUsage(); return; }

if (p.rank < 1)
{
PrintUsage(); return;
}

p.ProcessA();
}
else
{
PrintUsage();
}

}

public void ProcessA()
{

// Print out the preamble stuff
Console.WriteLine(
"\\documentclass{article}" + @"

% This is the output from LatexMatrixInverse 1.0 for a matrix with rank " + rank + @"

\usepackage{geometry}

% Note: you should probably use pdflatex to compiile this file.
% Other processors are known to have some issues with using
% 'geometry' to set paper size

% Adjust the page size here if output is wrapping in a bad way.
% Default is 8.5 x 11 in (Letter)
\geometry{papersize={8.5in,11in}}

%Import AMS Latex packages
\usepackage{amsmath, amssymb}
\setcounter{MaxMatrixCols}{" + (2 * rank + 1) + @"}

%Variable definition
");

// Define variables
//
// We have A = {a_i_j}, B = {b_j}
// Want to find solution to Ax=B by Gauss-Jordan elimination
//
// Variables format: X= row, Y = column
// \newcommand{ArXcY}{a_{X,Y}}
// \newcommand{BrX}{b_X}
// Have to do odd things with the numbers in the command names
// because latex doesn't support numbers in \newcommand

A = new string[rank, rank];
A_inverse = new string[rank, rank];
B = new string[rank];

StringBuilder A_definition = new StringBuilder("% Definition of initial A");
StringBuilder B_definition = new StringBuilder("% Definition of initial B\r\n");

for (int i = 0; i < rank; i++)
{
A_definition.AppendLine(@"
% A row " + (i + 1));

for (int j = 0; j < rank; j++)
{
A[i, j] = "\\AR" + NtoC(i + 1) + "C" + NtoC(j + 1);

A_definition.AppendLine("\\newcommand{\\AR" + NtoC(i + 1) + "C" + NtoC(j + 1)
+ "}{a_{" + (i + 1) + "," + (j + 1) + "}}");

A_inverse[i, j] = (i == j ? "1" : "0");

}

B[i] = "\\BR" + NtoC(i + 1);

B_definition.AppendLine("\\newcommand{\\BR" + NtoC(i + 1) + "}{b_{" + (i + 1) + "}}");
}

Console.WriteLine("\\begin{document}");

Console.WriteLine(A_definition);
Console.WriteLine(B_definition);
Console.WriteLine();
Console.WriteLine();

Console.WriteLine("LatexMatrixInverse 1.0 Output for rank " + rank + ", ShowIntermediateSteps is "
+ ShowIntermediateSteps.ToString() + ".");
Console.WriteLine();
Console.WriteLine();

Console.WriteLine("Given the following initial matrices:");

Console.WriteLine("\\begin{equation*}");
// Print out A
Console.Write("A = \\begin{pmatrix}");
for (int i = 0; i < rank; i++)
{
if (i > 0) Console.Write("\\\\");

for (int j = 0; j < rank; j++)
{
if (j > 0) Console.Write("&");

Console.WriteLine(A[i, j]);
}

}
Console.Write("\\end{pmatrix}");

// Print out B
Console.Write("B = \\begin{pmatrix}");
Console.Write(B[0]);
for (int i = 1; i < rank; i++)
{
Console.Write("\\\\ " + B[i]);
}
Console.Write("\\end{pmatrix}");
Console.WriteLine("\\end{equation*}");

Console.WriteLine(" We want to find $A^{-1}$ and $A^{-1} B$...");

StringBuilder array_fmt = new StringBuilder();

for (int j = 0; j < 2; j++)
{
for (int i = 0; i < rank; i++)
{
array_fmt.Append("c");
}

array_fmt.Append("|");
}
ArrayFormat = array_fmt.Append("c").ToString();

if (ShowIntermediateSteps)
{

Console.WriteLine("\\begin{equation*}");
Console.WriteLine(GetCombinedArrays());
Console.WriteLine("\\end{equation*}");

}

string multiplier = null;

//
// Build the upper triangle matrix in A
//
for (int i = 0; i < rank; i++)
{
multiplier = A[i, i];

B[i] = "\\frac{" + B[i] + "}{" + multiplier + "}";

if (A_inverse[i, i] != "0")
{
A_inverse[i, i] = "\\frac{" + A_inverse[i, i] + "}{" + multiplier + "}";
}

A[i, i] = "1";

//work with under the diagonal elements on the inverse

for (int j = 0; j < i; j++)
{
if (A_inverse[i, j] != "0")
{
A_inverse[i, j] = "\\frac{" + A_inverse[i, j] + "}{" + multiplier + "}";
}
}

//work with over the diagonal elements in the row for A

for (int j = i + 1; j < rank; j++)
{
A[i, j] = "\\frac{" + A[i, j] + "}{" + multiplier + "}";
}

if (ShowIntermediateSteps)
{
Console.WriteLine("\\begin{equation*}");
Console.WriteLine("\\xrightarrow{\\frac{1}{" + multiplier + "} R_{" + (i + 1) + "}}");
Console.WriteLine(GetCombinedArrays());
Console.WriteLine("\\end{equation*}");
}

for (int j = i + 1; j < rank; j++)
{
multiplier = A[j, i];

// Work with B
B[j] = B[j] + " - \\left (" + multiplier + " \\right ) \\left ( " + B[i] + "\\right )";

A[j, i] = "0";

for (int k = i + 1; k < rank; k++)
{
A[j, k] = A[j, k] + " - \\left (" + multiplier + " \\right ) \\left ( " + A[i, k] + "\\right )";
}

//Work with inverse matrix

for (int k = 0; k <= i; k++)
{
if (A_inverse[i, k] != "0")
{
A_inverse[j, k] = A_inverse[j, k] + " - \\left (" + multiplier + " \\right ) \\left ( " + A_inverse[i, k] + "\\right )";
}
}

if (ShowIntermediateSteps)
{
Console.WriteLine("\\begin{equation*}");
Console.WriteLine("\\xrightarrow{R_{" + (j + 1) + "} - \\left ( " + multiplier + "\\right ) R_{" + (i + 1) + "}}");
Console.WriteLine(GetCombinedArrays());
Console.WriteLine("\\end{equation*}");
}

}

}

//Now finish the elimination from the diagonal in A to the upper triangle of A
for (int i = rank - 1; i >= 0; i--)
{
for (int j = i - 1; j >= 0; j--)
{
//What's the multiplier?
multiplier = A[j, i];
A[j, i] = "0";

//Handle the B[j]
B[j] = B[j] + " - \\left (" + multiplier + " \\right ) \\linebreak \\left ( " + B[i] + "\\right )";

//Nothing to do with A in this case

//A inverse is going to be a full row operation
for (int k = 0; k < rank; k++)
{
A_inverse[j, k] = A_inverse[j, k] + " - \\left (" + multiplier + " \\right ) \\linebreak \\left ( " + A_inverse[i, k] + "\\right )";
}

if (ShowIntermediateSteps)
{

Console.WriteLine("\\begin{equation*}");
Console.WriteLine("\\xrightarrow{R_{" + (j + 1) + "} - \\left ( " + multiplier + "\\right ) R_{" + (i + 1) + "}}");
Console.WriteLine(GetCombinedArrays());
Console.WriteLine("\\end{equation*}");
}
}
}

if (FinalOutputAsMatrix)
{

Console.WriteLine("\\begin{equation*}");
// Print out A
Console.Write("A^{-1} = \\begin{pmatrix}");
for (int i = 0; i < rank; i++)
{
if (i > 0) Console.Write("\\\\");

for (int j = 0; j < rank; j++)
{
if (j > 0) Console.WriteLine("&");

Console.WriteLine(A_inverse[i, j]);
}

}
Console.Write("\\end{pmatrix}");

// Print out B
Console.Write("A^{-1}B = \\begin{pmatrix}");
Console.Write(B[0]);
for (int i = 1; i < rank; i++)
{
Console.Write("\\\\ " + B[i]);
}
Console.Write("\\end{pmatrix}");
Console.WriteLine("\\end{equation*}");

}
else
{
Console.WriteLine();
Console.WriteLine();
Console.WriteLine("The solution for $A^{-1}$ follows...");

Console.WriteLine();

for (int i = 0; i < rank; i++)
{
for (int j = 0; j < rank; j++)
{
Console.WriteLine("Row " + i + ", Column " + j + " of $A^{-1}$ is");
Console.WriteLine("\\begin{equation*}");
Console.WriteLine(A_inverse[i, j]);
Console.WriteLine("\\end{equation*}");

}
}

for (int i = 0; i < rank; i++)
{
Console.WriteLine("Row " + i + " of $A^{-1} B$ is");
Console.WriteLine("\\begin{equation*}");
Console.WriteLine(B[i]);
Console.WriteLine("\\end{equation*}");

}

}

//End of file
Console.WriteLine("\\end{document}");

}

public string GetCombinedArrays()
{
StringBuilder output = new StringBuilder();

output.Append("\\left ( \\begin{array}{" + ArrayFormat + "}");

for (int i = 0; i < rank; i++)
{
if (i > 0) { output.AppendLine("\\\\"); }

for (int j = 0; j < rank; j++)
{
if (j > 0) output.Append("&");
output.AppendLine(A[i, j]);
}

for (int j = 0; j < rank; j++)
{
output.Append("&");
output.AppendLine(A_inverse[i, j]);
}

output.Append("&");
output.Append(B[i]);

}

output.Append("\\end{array} \\right )");
return output.ToString();
}

public string NtoC(int input)
{
return NtoC(input.ToString());
}

public string NtoC(string input)
{
return input.Replace("0", "a")
.Replace("1", "b")
.Replace("2", "c")
.Replace("3", "d")
.Replace("4", "e")
.Replace("5", "f")
.Replace("6", "g")
.Replace("7", "h")
.Replace("8", "i")
.Replace("9", "j");

}

public static void PrintUsage()
{
Console.WriteLine(@"
LatexMatrixInverse C#.Net Edition 1.0
Developed by Mike Burr
This application is provided without any warranties, express or implied.

Usage: LatexMatrixInverse <rank>

<rank>: The rank of the matrix. Also the number of the diagonal elements in a square matrix.
ex. A 2x2 matrix has rank 2, a 3x3 matrix has rank 3, etc...

The latex code for the matrix inverse and corresponding solution vector
of a matrix with given rank will be written to standard out along with supporting work
code (ex. row transformations and additions). The end user should then adjust the matrix
values in the resulting latex code to get the desired result for their specific problem.

The matrix inversion is performed using Gauss-Jordan elimination.

Note that singular matrices still won't be invertible, but matrices with inverses
should be invertible using the latex code generated with the occasional initial
row substitution. Simplification of the resulting formulas are the end user's
responsibility.

A recommended use of the application would be to redirect the output to a file.
ex. The inverse for a 3x3 matrix written to m3x3.tex:

LatexMatrixInverse 3 > m3x3.tex
");
}
}
}

-->