7
$\begingroup$

I am a regular user of Mathematica, Julia, and MATLAB but I am looking for something different. The problem I am trying to solve in Mathematica only requires (dense) linear algebra to specify but is of a size such that Mathematica takes too long (and seems to automatically halt the computation after like 18 hours, when it should finish in a week or so). I am looking for a symbolic solution and not a numeric solution (to then use a machine learning approach on the resulting equations). I was wondering if anyone knew if there are more efficient and parallel symbolic packages. Checking around it looks like Singularity or Machaly2 may be what's right for this problem?

Edit: As requested here is more detail. It's as straight-forward as it sounds. I need to solve for the symbolic solution of a system of 4x4 matrices and 1x4 vectors. I am able to succinctly write down the exact solution (attached as a photo, hard to format it right in the editor...),

Exact Solution

but as you can see it's a little unwieldy (then it gets squared and simplified). I want to use the resulting polynomial (in terms of the matrix coefficients) in numerical computations, but I don't want to have to solve this system of equations each time the matrices change since this would increase the computational cost immensely.

So at a higher level it's really simple, and coding Mathematica to get it to try this computation is simple, but it is taking a lot of computing time (largely because Mathematica's engine is single threaded). In the simple case when the $A$ and $B$ matrices are lower triangular, I Mathematica was able to solve it after about 10 minutes. But getting rid of this assumption on even one matrix increases the complexity by a lot.

Let me lastly note that I have computing power at my disposal. I do high-performance computing, and so I have a few nodes available for a few hundred CPU cores and about a terabyte of RAM to let it churn for a month (I was hoping to change each matrix from lower triangular to full one by one to see how far I can get). However, I don't know what program/language will do this kind of "HPC Symbolic Computing". I was hoping there was some computational algebra package that has the tools to do this.

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ If you described your problem, perhaps something useful could be said, otherwise, your question is too general. $\endgroup$
    – Igor Rivin
    Jan 26, 2016 at 23:44
  • 1
    $\begingroup$ What do you mean symbolic and linear algebra? So, you want to do gaussian elimination, but with some unknown entries? $\endgroup$
    – Per Alexandersson
    Jan 27, 2016 at 0:00
  • 1
    $\begingroup$ Edited with more detail. I mean symbolic computation in the standard sense, and I mean linear algebra as in the equation to simplify can be expressed in terms of linear algebra objects and operations. All of the matrix coefficients are unknowns, and am looking for a solution in terms of the coefficients. $\endgroup$
    – Chris Rackauckas
    Jan 27, 2016 at 0:15
  • $\begingroup$ Did you try to ask on the dedicated Mathematica SE site? $\endgroup$
    – kjetil b halvorsen
    Jan 27, 2016 at 13:02
  • 1
    $\begingroup$ I assumed this would be off-topic on the Mathematica SE site since it's not about Mathematica. Due to the Mathematica's general overhead and not having a multithreaded kernal, I thought it would make more sense to ask other mathematicians who encountered the same problem and had switch from Mathematica. I know some reasearchers in things like computational group theory have specialized packages for solving computationally hard problems (GAP), so I was hoping to get advice from some of those folks here. However, it seems I may need to seek them out on the more computational forums. $\endgroup$
    – Chris Rackauckas
    Jan 27, 2016 at 17:08

2 Answers 2

Reset to default
4
$\begingroup$

I found out that the best solution is actually SymEngine.jl in Julia. Using Julia as the driver ends up increasing the speed and reducing the memory requirement, making it more efficient than using SymEngine in Python (or using SymPy). This ended up being much more efficient than Mathematica, and the few things I experimented with in Sage. Also, Julia's operator overloading lets you write a "generic" matrix function and then it will automatically compile a function which gets very good performance on symbolic expressions, and this makes it easy to extend the symbolic libraries as necessary without losing performance. This ended up being a very good solution to brute force some symbolic math to setup further analysis.

Improve this answer
$\endgroup$
2
$\begingroup$

If $M$ is an invertible matrix and $\bf v$ is a vector, Cramer's rule expresses $M^{-1} \bf v$ as $\det(M)^{-1} \bf w$ where $\det(M)$ and the entries of $\bf w$ are polynomials in the entries of $M$ and $\bf v$. In general, you can't make it any simpler than that. The more complicated the expression, the less likely you are to have cancellations.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I am not expecting many cancellations. In the case where I use lower triangular matrices, I get something which simplifies to about 10,000 terms. I expect to have a lot more if I compute the full equation out. $\endgroup$
    – Chris Rackauckas
    Jan 27, 2016 at 0:28
  • 2
    $\begingroup$ Yes, expanding everything out is likely to produce a combinatorial explosion. Definitely not recommended. $\endgroup$
    – Robert Israel
    Jan 27, 2016 at 0:30
  • 1
    $\begingroup$ I was hoping to solve it in an HPC manner, like using a few hundred CPUs and a tarabyte of memory and coming back a month later. For numeric computations there are many ways to do this. For symbolic I don't know of a package/language which does this, and I wouldn't want to try from scratch. $\endgroup$
    – Chris Rackauckas
    Jan 27, 2016 at 0:34
  • 1
    $\begingroup$ As I said, the symbolic solution of a $4 \times 4$ invertible system is expressible using Cramer's rule. Each $4 \times 4$ determinant is the sum of $24$ terms, where each term is $\pm 1$ times a product of four factors. In that form, the expression is very easy to evaluate. No need for hundreds of CPU's or terabytes of memory. If you try to combine a bunch of these and expand everything out, you get a nightmare. So don't expand everything out. $\endgroup$
    – Robert Israel
    Jan 27, 2016 at 20:53
  • 1
    $\begingroup$ Thanks for the clarification. That might work. It's still less than ideal since every step of the machine learning algorithm I will have to recompute anything where the $A$'s and $B$'s are involved (for multiple $\mu$ and $\sigma$), but your method definitely puts it in the "computable" range. I threw the Mathematica code on the HPC to just churn on one node and see if that finishes before I get it solved with your method. However, I am going to leave the topic open to see if anyone answers with some package they are working on since that would be useful for similar problems in the future. $\endgroup$
    – Chris Rackauckas
    Jan 28, 2016 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.