# Compiling LinearSolve[] or creating a compilable procedural version of it

Earlier today I had a discussion with a representative at Premier Support about the 2 questions I've asked here over the past couple of days:

Neither the conversation nor the answers to the above questions have produced a solution as good as compiling my function to C code. This leaves me with some problems as the function I need to deploy uses 4 high-level Mathematica functions which don't (or maybe don't readily) compile:

CholeskyDecomposition[]
IdentityMatrix[]
LinearSolve[]
Simplify[]


Yesterday I made some progress developing/adapting solutions for procedural versions of CholeskyDecomposition[] and IdentityMatrix[], so these at least should compile.

Given the anticipated use of my function, maybe I can get away without using Simplify[].

This would leave me still needing a procedural equivalent to LinearSolve[] or some means of compiling it.

Maybe simpler because, I don't need a myLinearSolve[] to do everything Mathematica's version does.

Interestingly, in my conversation with the guy from support, he suggested that given the limited nature of what I wanted to do, he thought that setting options on LinearSolve[] should give me a reduced scope version of LinearSolve[] that Mathematica could compile. He couldn't specify which options.

So let's explore this.

A typical use of LinearSolve[] in my current function looks like this:

c = {0.516344, 0.287671, 0.216344, 0.176796, 1};
A = {{0, 1, 1/2, 1/2, 1/2}, {0, 0, Sqrt[3]/2, 1/(2 Sqrt[3]), 1/(2 Sqrt[3])}, {0, 0, 0, Sqrt[2/3], 1/(2 Sqrt[6])}, {0, 0, 0, 0, Sqrt[5/2]/2}, {1, 1, 1, 1, 1}};

LinearSolve[A, c]

{0.173339, 0.206027, 0.187944, 0.209058, 0.223631}


No symbolic processing. No imaginary numbers. No SparseArray objects. A always a square matrix; c always a vector; and its usage always outputs a real vector.

Question 1 -- Does anyone have any insight into how to set options for the above use of LinearSolve[] so it will compile or if this is even possible?

Note: Oleksandr's answer in Why is MainEvaluate being used when LinearSolve can be compiled? may have some bearing on this.

...

If no way forward comes from the above, I may still have a chance to implement a limited procedurally based proceduralLinearSolve[].

If very distant memory serves me, one can use the inverse of a square matrix for some problems addressable with LinearSolve[]. This might give me a more specific way forward except Mathematica's implementation of it, Inverse[], doesn't fall on the list of compilable functions either.

I have found some code for calculating the inverse of a square matrix. Not even certain of the language, but I can probably follow its logic and port it to a proceduralInverse[] function in Mathematica.

This background takes me to my second question...

Question 2 -- Does creating a proceduralLinearSolve[] or proceduralInverse[] seem like a viable way to go (has anyone tried this kind of thing and succeeded) and/or can anyone point me in the direction of resources or solutions that could help me do this?

-
These functions are not compilable because they are calls to the Intel MKL. Don't waste your time coding them youself, and for a small function there is no reason to use a high-performance BLAS library either. Use LAPACK/CLAPACK (example) or Eigen. The response of the support guy is curious and, I think, incorrect. – Oleksandr R. Jan 14 '13 at 20:27
Depending on how systematic the elements you are passing to LinearSolve, you could just solve it once and compile the resulting formula for the solution instead of trying to compile the function itself. – jVincent Jan 14 '13 at 21:51
Using Inverse is not advisable. Typically, one uses LinearSolve to avoid inverting matrices. – asim Jan 14 '13 at 23:06
@OleksandrR. -- Thanks for the suggestions. I have hoped that I could rely on Mathematica for all my development needs. It just seems odd that with all the amazing things Mathematica does, it gets tripped up on deploying solutions. – Jagra Jan 14 '13 at 23:58
@Jagra You shouldn't be concerned with the return value. The question is what your input matrix looks like, for instance if it was {{a,b},{0,d} with rhs {r1,r2} where you know that a and d are non-zero, then you can call it once to see that your solution will always be: {(d r1 - b r2)/(a d), r2/d}. Now this method of cause will not work if you have an arbitrary matrix, but it you have a pattern and only some inputs, you can just solve once and compile the resulting definition. – jVincent Jan 15 '13 at 0:14

In many cases you don't really need to have LinearSolve compiles, since you are dealing for instance with only a specific type of matrices, or maybe just want to solve the same system with different right hand sides. In such cases you can evaluate the definiions once and compile the resulting expressions. Here is a rough example where I just injected two unkowns into your matrix.

c = {0.516344, 0.287671, 0.216344, 0.176796, 1};
A = {{0, b, 1/2, 1/2, 1/2}, {0, 0, Sqrt[3]/2, a/(2 Sqrt[3]),
1/(2 Sqrt[3])}, {0, 0, 0, Sqrt[2/3], 1/(2 Sqrt[6])}, {0, 0, 0, 0,
Sqrt[5/2]/2}, {1, 1, 1, 1, 1}};

fun[b_, a_] = LinearSolve[A, c]
funC = Compile[{{a, _Real}, {b, _Real}}, Evaluate[fun[a, b]]]

CompiledFunction[{a,b},
Block[{Compile$2,Compile$9},
Compile$2=1/a; Compile$9=1. b;
{0.0696861 Compile$2 (-2.4565+4.44393 a-0.5 b+1. a b), 0.0348431 Compile$2 (4.913 +Compile$9), -0.0696861 (-3.69701+Compile$9),
0.209058,
0.223631}] ,-CompiledCode-]


Now however when doing this it's important to remember that you need to check your inputs yourself, for instance in the case of A={{p1,0},{0,p2}}, you would need to check that neither p1 or p2 are zero.

-

Rob Knapp has some excellent publicly available notebooks for some of the things you might want. Here are three in particular that may be useful.

http://www.mathematica-journal.com/issue/v7i4/features/knapp/

http://library.wolfram.com/infocenter/Conferences/288/

http://library.wolfram.com/infocenter/Conferences/7968/

Among other things they give explicit implementations of Cholesky decomposition and linear solving, very likely suitable to run through Compile with C as target.

-
Thanks or the links. I've started to breath easier. – Jagra Jan 14 '13 at 22:05
Daniel -- I found Rob's Cholesky Decomposition and it looks much more robust than what I've come up with so far. He runs through a number of solving methods, did you have an idea of which of the one's he demonstrates would apply? – Jagra Jan 15 '13 at 1:04
Not offhand. Probably will depend on details of your intended application, typical inputs, etc. – Daniel Lichtblau Jan 15 '13 at 14:26

Using as basis the great resource for core numerical algorithms below, I managed to implement a compiled linear solve which doesn't call MainEvaluate (so quite fast).

I needed a linear solve for an optimization where the objective function requires inverting matrices, I was hesitating to use C++, but I preferred to stay in Mathematica.

Resources

Example

A=RandomReal[UniformDistribution[],{3,3}];
b={1,2,3};
A.LinearSolveCompiled[A,b]


Implementation

LUCompiled =
Compile[{{A0,_Real,2}},
Module[{i,j,k,p,r,A,n,L,c,U, P},

A=A0;
n=Length@A;
r=Table[j,{j,n}];

For[p=1,p<=n-1,p++,
For[j=p+1,j<=n,j++,
If[Abs[A[[r[[j]],p]]]>Abs[A[[r[[p]],p]]],
r[[{j,p}]]=r[[{p,j}]];
];
];

For[k=p+1,k<=n,k++,
A[[r[[k]],p]]=A[[r[[k]],p]]/A[[r[[p]],p]];

For[c=p+1,c<=n,c++,
A[[r[[k]],c]]=A[[r[[k]],c]]-A[[r[[k]],p]]*A[[r[[p]],c]];
];
];
];

L=Table[0.,{n},{n}];
For[i=1,i<=n,i++,
L[[i,i]]=1.;
];

P = L;
P = P[[r]];

U=A[[r]];

For[i=1,i<=n,i++,
For[j=1,j<=i-1,j++,
L[[i,j]]=A[[r[[i]],j]];
U[[i,j]]=0.;
];
];

Join[L,U,P]
]
];

BackSubCompiled=
Compile[{{a,_Real,2},{b,_Real,1}},
Module[{i,j,x,sum,n},

n=Length@a;
x=Table[0.,{n}];

For[i=n,1<=i,i--,

sum=0.;
For[j=i+1,j<=n,j++,
sum+=a[[i,j]]*x[[j]];
];

x[[i]]=(b[[i]]-sum)/a[[i,i]];
];

x
]
];

ForeSubCompiled=
Compile[{{a,_Real,2},{b,_Real,1}},
Module[{i,j,x,sum,n},

n=Length@a;
x=Table[0.,{n}];

For[i = 1, i <= n, i++,

sum=0.;
For[j=1,j<=i-1,j++,
sum+=a[[i,j]]*x[[j]];
];

x[[i]] = (b[[i]] - sum)/a[[i,i]]
];

x
]
];

(*Solves x so that A.x = b*)
LinearSolveCompiled=
Compile[{{A,_Real,2},{b,_Real,1}},
Module[{n,luResult, L, U, P, x1, x2},

n=Length@A;
luResult=LUCompiled[A];

L=luResult[[;;n]];
U=luResult[[n+1;;2n]];
P=luResult[[2n+1;;3n]];

x1=ForeSubCompiled[L,P.b];
x2=BackSubCompiled[U,x1];

x2
]
,
CompilationOptions->{"InlineCompiledFunctions"->True,"InlineExternalDefinitions"->True}
];


Checking for calls to MainEvaluate

FastCompiledFunctionQ[function_CompiledFunction]:=
(
Needs["CompiledFunctionTools"];
StringFreeQ[CompiledFunctionToolsCompilePrint@function, "MainEvaluate"]
);

LinearSolveCompiled // FastCompiledFunctionQ (*True*)

-