1
$\begingroup$

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

Improve this question
$\endgroup$
5
  • $\begingroup$ You should examine what NullSpace of L is. See e.g. closely related post: All possible solutions to the Matrix Equation (free variables appearing). $\endgroup$
    – Artes
    Commented Jul 9, 2014 at 13:09
  • $\begingroup$ 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. $\endgroup$
    – dreichler
    Commented Jul 9, 2014 at 13:27
  • $\begingroup$ 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. $\endgroup$
    – Artes
    Commented Jul 9, 2014 at 13:35
  • $\begingroup$ 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? $\endgroup$
    – dreichler
    Commented Jul 9, 2014 at 16:00
  • $\begingroup$ The answer given points out one of possible issues, this seems quite reasonable since no further details of the problem were exposed. $\endgroup$
    – Artes
    Commented Jul 9, 2014 at 16:20

1 Answer 1

Reset to default
5
$\begingroup$

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

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – dreichler
    Commented Jul 9, 2014 at 15:55
  • $\begingroup$ 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. $\endgroup$
    – Dr. Wolfgang Hintze
    Commented Jul 9, 2014 at 17:30
  • $\begingroup$ 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. $\endgroup$
    – dreichler
    Commented Jul 10, 2014 at 8:12
  • $\begingroup$ 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. $\endgroup$
    – Dr. Wolfgang Hintze
    Commented Jul 10, 2014 at 12:21
  • $\begingroup$ 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. $\endgroup$
    – dreichler
    Commented Jul 11, 2014 at 10:52

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.