Skip to main content

matrix - 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?



Answer



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.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more