# LinearSolve failing to find solution to system of linear equations of 64 variables

I've got a 64x64 Matrix L for which I want to find the solution to the matrix equation L . x == rho.

L is defined as L[a, b, c, d, e, f] with assumptions Element[{a, b, c, d, e, f}, Reals].

As long as I put specific numbers in L, e.g. L[1, 1, 1, 1, 1, 1], LinearSolve[L, rho] (with rho something like rho = {1, 0, 0, ..., 0}) finds the solution within a fraction of a second. But when making even one argument symbolic, say L[a, 1, 1, 1, 1, 1], LinearSolve finds no solution (with rho as defined above):

LinearSolve::nosol: Linear equation encountered that has no solution.

But there is a solution, at least all values of a I tried manually. And a is only added linearly in only some of the components of L. And even Det[L] is nonzero, so there should be a solution. I assume I'm doing something pretty wrong.

Can anyone help me?

Edit:

Here is a sample matrix L[a,b,c,d,e,f] with fewer dimensions(5x5) where LinearSolve works(probably not the best sample because the big L matrix is more or less spare and this example doesn't look so):

L[a,b,c,d,e,f]={{1, 0, 0, 0, 0}, {0, (0. + 2. I) c, -((I dSin[f])/Sqrt[3]), -((I d Cos[f])/Sqrt[3]), 0},
{(I d Cos[f])/Sqrt[3], -((I d Sin[f])/Sqrt[3]), -11.565 - (0. + 1. I) a +(0. + 0.666667 I) c, 0, -(1/2) I e Sin[f]},
{-((I d Sin[f])/Sqrt[3]), -((I d Cos[f])/Sqrt[3]), 0, -11.565 - (0. + 1. I) a + (0. + 1.33333 I) c, 0},
{0, 0, -(1/2) I e Sin[f], 0, -0.2 - (0. + 1. I) a + (0. + 1. I) b - (0. + 0.2 I) c}}


i is no variable but the imaginary unit

• You should examine what NullSpace of L is. See e.g. closely related post: All possible solutions to the Matrix Equation (free variables appearing). Jul 9, 2014 at 13:09
• NullSpace[L] with one variable in L, so L[a,1,1,1,1,1] gives a very very long expression with exponents up to a^22 and more. Sorry, I didn't understand how Nullspace[] brings me any closer to the solution of the matrix equation. I tried to understand what's going on in the thread you linked but didn't fully understand. Jul 9, 2014 at 13:27
• Understanding linear algebra is indispensable to proceed further. I suggest reproducing the issue you encounter with a smaller example, otherwise I doubt anyone will help you. Jul 9, 2014 at 13:35
• I tried to reproduce the error with smaller matrices but it didn'nt happen(i did not try further than 5x5). LinearSolve evaluates a symbolic solution as expected. I don't get what's the principal difference between a 5x5 and a 64x64 matrix. Why is it possible in a smaller one? Jul 9, 2014 at 16:00
• The answer given points out one of possible issues, this seems quite reasonable since no further details of the problem were exposed. Jul 9, 2014 at 16:20

The problem here can arise because of numerical underflow which appears for sufficiently large dimension of the problem. Some numerically very small number multiplies the parameter "a" and therefore "a" does not appear in the "solution".

Consider a simple example

Define the matrix m (fill it with random numbers, here exponentially distributed)

In[263]:= n = 10;
m = Table[-Log[Random[]], {i, 1, n}, {j, 1, n}];


Now replace one element by the parameter "a":

In[265]:= m[[1, 1]] = a

Out[265]= a


Check the determinant

In[266]:= Det[m] // Distribute

Out[266]= -364.727 + 119.62 a


Similarly define the right hand side

In[267]:= y = Array[-Log[Random[]] &, n]

Out[267]= {0.335897, 0.0365035, 1.09626, 0.0109352, 2.37241, 0.364962, 0.682719, \
1.67369, 1.05203, 1.07763}


Solve the equation

In[269]:= ls = FullSimplify[LinearSolve[m, y]]

Out[269]= {2.85656*10^-15 - 5.42695*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a),
1.47799 + 2.29225*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a), -1.48418 +
2.10951*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a),
0.00418502 + 3.91494*10^-169/(
6.28005*10^-169 - 2.05968*10^-169 a), -3.52696 + 5.26631*10^-169/(
6.28005*10^-169 - 2.05968*10^-169 a), -0.0479913 - 3.29329*10^-169/(
6.28005*10^-169 - 2.05968*10^-169 a), -1.58069 + 2.58972*10^-169/(
6.28005*10^-169 - 2.05968*10^-169 a),
0.0154976 + 3.84753*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a),
1.00619 - 3.66389*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a),
1.8942 - 5.07079*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a)}


A typical term in which "a" appears is

5.42695*10^-169/(6.28005*10^-169 - 2.05968*10^-169 a)


We can see that there appears a very small factor 2.05968*10^-169 before "a".

But this is of course spurious because it can be reduced by dividing by the numerator. Let's do this using a pure function

f = 1/Distribute[Denominator[#]/Numerator[#]] &


Which gives something looking innocently:

In[270]:= f[5.426948941128063*^-169/(
6.280047143587323*^-169 - 2.059678646225578*^-169 a)]

Out[270]= 1/(1.1572 - 0.379528 a)


We can extend this to the whole list ls:

In[271]:= ls1 = Plus @@ (f /@ List @@ #) & /@ ls

Out[271]= {2.85656*10^-15 + 1/(-1.1572 + 0.379528 a),
1.47799 + 1/(2.73969 - 0.898542 a), -1.48418 + 1/(2.97702 - 0.976379 a),
0.00418502 + 1/(1.60413 - 0.526108 a), -3.52696 + 1/(
1.1925 - 0.391105 a), -0.0479913 + 1/(-1.90692 + 0.625417 a), -1.58069 + 1/(
2.42499 - 0.795328 a), 0.0154976 + 1/(1.63223 - 0.535325 a),
1.00619 + 1/(-1.71404 + 0.562157 a), 1.8942 + 1/(-1.23847 + 0.406185 a)}


Ok, done. Now you can put in specific values vor "a":

In[222]:= ls1 /. a -> 10

Out[222]= {-0.114464, 0.335937, -0.13769, -0.0448387, 0.133666, 0.207523, 0.05006, \
-0.427603, 0.0156447, 0.198275, -0.0961547, -0.057508, 0.191488, 0.336546, \
-0.359727, 0.301349}

In[153]:= ls /. a -> 1.

Out[153]= {16.4706, 1.92368, -5.40107, 0.527357, 3.18836, 1.39021, 0.36902, 1.17988, \
-3.67465, -0.996117, -5.23853, -2.44246, -1.243, -0.426711, -0.939764, 0.}

In[154]:= ls /. a -> -1

Out[154]= {16.4706, 1.92368, -5.40107, 0.527357, 3.18836, 1.39021, 0.36902, 1.17988, \
-3.67465, -0.996117, -5.23853, -2.44246, -1.243, -0.426711, -0.939764, 0.}
`

I suggest playing around with this model, varying the dimension. You will notice that from n~=15 on there will be no parameter "a" any more in the solution because of numerical underflow.

Best regards, Wolfgang

• Hey, very interesting. I don't know whether this is the problem in my case because when I give a the value a=0(so make it vanish) the LinearSolve[L[0,1,1,1,1,1,1],rho] does his job and gives me a solution. Jul 9, 2014 at 15:55
• Please provide an example of your data (in Code Format). Let the dimension be much smaller than 64. It will be my pleasure to investigate what's going on. Jul 9, 2014 at 17:30
• Thanks for the help so far. I edited my original post. Now there is some example of my data. But LinearSolve works with the example. If you're interested in the full matrix/mathematica code just let me know and I will send it to you. Jul 10, 2014 at 8:12
• We can take up the topic on Monday. As I said before and have shown in my model the problem will appear for suffiently high dimension (see my example). You can try your examples on your own, increase the dimension and find the "critical" one. I guess it will be Close to 15. Jul 10, 2014 at 12:21
• The critical one is 60. Everything works just fine incl. 59. Above that border mathematica claims there is no solution, but there is, i tried it at least for 1000 values of a in the interval -100,100. (i'm not interested in other ranges anyway). I would be very very glad if you could help me investigate any further. But now i'm in holiday for two weeks and cannot check the thread and your answers. But if you write me something i will get it in two weeks, thread is bookmarked. Jul 11, 2014 at 10:52