# Solving determinant of a Kronecker product of matrices gives a numerical error - why?

I am doing the following steps (code at the end of the post):

1. I start with a 2x2 matrix (smatrix), which is a function of a single variable (u2). I want to set the determinant of this matrix (smatrix) to zero, therefore solving for u2.

2. There are two real solutions (sols), which Solve (or NSolve) delivers. These are exact, as the determinant is the order 2 polynomial in the variable (u2). It might be important that the coefficients of this polynomial are on the order of 10^-10.

3. I make a 4x4 matrix (matrix2) which is just the earlier 2x2 matrix two times on the diagonal (KroeneckerProduct). Its determinant is then the determinant of the first matrix squared.

4. I solve the new determinant to be zero, now getting different solutions. They are complex. Moreover, when plugging in the old and the new solutions into the new determinant, I see that the new solutions apparently are more accurate (by 10^4!).

I have tried playing around with global WorkingPrecision, but it does not seem to change anything. Can somebody explain to me what goes wrong here? Moreover, what should I change in my code to retain the exact results?

Thanks a lot!

smatrix = {{1. - 2.96392/u2,  0.0000196744/ u2},
{1. - 2.96392/u2, -2.11737*10^-10 + 0.0000196746/u2}};
sols = (smatrix // Det) // Solve[# == 0, u2] &


{{u2 -> 0.987752}, {u2 -> 2.96981}}

(smatrix // Det) /. sols[[1]]


0.

matrix2 = KroneckerProduct[({
{1, 0},
{0, 1}
}), smatrix];
(matrix2 // Det) // Solve[# == 0, u2] &


{{u2 -> 0.987752 - 4.12625*10^-6 I} , {u2 -> 0.987752 + 4.12625*10^-6 I} , {u2 -> 2.96981 - 0.000012998 I} , {u2 -> 2.96981 + 0.000012998 I} }

(matrix2 // Det) /. %[[1]]


3.85186*10^-34 - 1.84553*10^-36 I

(matrix2 // Det) /. sols[[1]]


3.15024*10^-30

-
You might want to give people your polynomial so they can replicate your issue. Also, you say your solutions for u2 are "exact" but clearly they are Reals, even though they could be "exact" in some other sense. If you use decimal notation (i.e. Reals) in your smatrix or your unspecified polynomial, you'll get inexact/numerical results from then on. – Kellen Myers Aug 26 '14 at 17:17
What specifically is smatrix? I've no idea how to replicate this without having the specific matrix. – Daniel Lichtblau Aug 26 '14 at 17:44
Thanks for your comments! Sorry for not including smatrix to start with, here it goes: In[288]:= smatrix Out[288]= {{1 - 2.96392/u2, 0.0000196744/ u2}, {1. - 2.96392/u2, -2.11737*10^-10 + 0.0000196746/u2}} – ffc Aug 26 '14 at 19:08
In case it gets buried at the bottom of my answer, I'll reiterate here that you really ought to be using Chop in these cases. – Kellen Myers Aug 26 '14 at 20:40
I thank all of you guys: Daniel Lichtblau, hieron, and Kellen Myers! I really appreciate the answers. The situation is much clearer for me now. I don't know how upvoting and accepting works, but I will try to do as much good as possible:) – ffc Aug 29 '14 at 22:55

This is really totally normal. Your characteristic polynomial has coefficients ~10^(-10) and your solutions sometimes seem to evaluate to f(sol)= 10^(-20) instead of 0. The precision with which Mathematica will NSolve an equation may be different than the care it takes when you plug in various floating-point reals back into the function they might satisfy. For example:

 ReplaceAll[z^100, NSolve[z^100 == 2][[55]]]


will give you 2. - 6.72179*10^-15 I, which looks wrong... but we all know it's probably fine, and the 55th solution to this equation is not somehow buggy or broken.

If you do the following:

Clear[smatrix]
smatrix[u2_] = {{1 - 2.96392/u2,
0.0000196744/u2}, {1. - 2.96392/u2, -2.11737*10^-10 +
0.0000196746/u2}}


And define:

f[u2_] = Expand[u2^2 Det[smatrix[u2]]]


I have multiplied by u2^2 to make this into a polynomial, so that our understanding of size of coefficients" makes more sense. If you don't do this, you'll get even worse answers. (See below.)

You can examine three sets of solutions:

sol2 = Solve[f[u2] == 0]
sol4 = Solve[f[u2]^2 == 0]
sol4a = Solve[u2^4 Det[KroneckerProduct[smatrix[u2], IdentityMatrix[2]]] == 0, u2]


Notice that sol4 and sol4a are supposed to be the same as sol2 (but with double multiplicity). On the other hand, sol4a gives us this error:

Solve::ratnz: Solve was unable to solve the system with inexact coefficients.
The answer was obtained by solving a corresponding exact system and numericizing
the result.


This is what user hieron suggests, to simply declare "These floating point real numbers are exact" but that won't actually lead you to an exact answer.. you're still abusing numbers that are not exact -- and never will be, without changing whatever it is you did to obtain smatrix in the first place. You can lie to Mathematica by rationalizing them, but that's just a dodge around this error, not a real fix. Dare I say, that is the wrong thing to do.

You can plug in:

Det[smatrix[u2]] /. sol2
Det[smatrix[u2]] /. sol4
Det[smatrix[u2]] /. sol4a


The first two cases give me all zeros (i.e. 0. the floating point real number zero). The third gives me very small complex numbers, indicating that there is some numerical error here. It's even worse if I try:

Map[Det, smatrix[u2] /. sol2]
Map[Det, smatrix[u2] /. sol4]
Map[Det, smatrix[u2] /. sol4a]


Doing the order-of-operations differently, the determinant coming after plugging in the solutions, it's even worse. But still, the error is 10 orders of magnitude smaller than the (already quite small) coefficients of my original f[u2] that I was solving.

Changing this all slightly, you can get very different results -- less accurate -- by simply failing to convert your functions into polynomials before solving, again indicating that this "inaccuracy" stems from numerical solving and, more importantly, that it is truly negligible:

g[u2_] = Expand[Det[smatrix[u2]]]
s2 = Solve[g[u2] == 0]
s4 = Solve[g[u2]^2 == 0]
s4a = Solve[
Det[KroneckerProduct[smatrix[u2], IdentityMatrix[2]]] == 0, u2]
Det[smatrix[u2]] /. s2
Det[smatrix[u2]] /. s4
Det[smatrix[u2]] /. s4a
Map[Det, smatrix[u2] /. s2]
Map[Det, smatrix[u2] /. s4]
Map[Det, smatrix[u2] /. s4a]


All of these errors (the last six lines, which in theory should all be zero) will be larger than before, but still relatively small (many orders of magnitude smaller than the coefficients). In fact, you're off by 2-4 orders of magnitude, since that is the power of u2 by which we failed to multiply. It may even make sense, depending on your numbers, to also multiply through by some power of 10. For example, if you have 40 digits of precision but your constants are ~10^(-10), why not multiply the equation by u2^2 and also 10^10? This may help you (and Mathematica) differentiate between which numbers are important and which are really error terms.

The use of a polynomial vs. a rational function vs. a function defined by a matrix vs. the square of some function vs. the function defined by a matrix product.... all of these are mathematically equivalent, but in each case we can make poor choices that increase the numerical error slightly. However, this error is normal.

Notice what happens when we do:

Chop[Det[smatrix[u2]] /. sol2]
Chop[Det[smatrix[u2]] /. sol4]
Chop[Det[smatrix[u2]] /. sol4a]
Chop[Map[Det, smatrix[u2] /. sol2]]
Chop[Map[Det, smatrix[u2] /. sol4]]
Chop[Map[Det, smatrix[u2] /. sol4a]]
Chop[Det[smatrix[u2]] /. s2]
Chop[Det[smatrix[u2]] /. s4]
Chop[Det[smatrix[u2]] /. s4a]
Chop[Map[Det, smatrix[u2] /. s2]]
Chop[Map[Det, smatrix[u2] /. s4]]
Chop[Map[Det, smatrix[u2] /. s4a]]


Zeros across the board. That's the purpose of Chop. Use it judiciously and carefully though, in particular along with the recommendation to normalize your coefficients to be ~1 instead of ~10^(-10).

-

I am not getting the same results in version 10. Also I see no problem with the results already indicated, and I'll say a bit about that as we proceed.

Here is what I obtain in version 10.

smatrix = {{1 - 2.96392/u2,
0.0000196744/u2}, {1. - 2.96392/u2, -2.11737*10^-10 +
0.0000196746/u2}};
det = Det[smatrix];
sols = Solve[det == 0, u2]

During evaluation of In[56]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(* Out[58]= {{u2 -> 0.944568025431}, {u2 -> 2.96391999999}} *)


Using instead NSolve:

nsols = NSolve[det == 0, u2]

(* Out[70]= {{u2 -> 2.96391999999}, {u2 -> 0.944568025431}} *)


Compute residuals:

det /. nsols

(* Out[77]= {5.16987882846*10^-26, 1.03397576569*10^-25} *)


Is this a viable result? The answer is that it is, but this residual check is not in and of itself conclusive; it does not take into account scaling of either inputs or results.

Now repeat for the Kronecker product matrix.

matrix2 = KroneckerProduct[({{1, 0}, {0, 1}}), smatrix];
det2 = Det[matrix2]

(* Out[80]= (1/(u2^4))1. (3.5139287066*10^-19 - 9.8114232253*10^-19 u2 +
9.35903247297*10^-19 u2^2 - 3.50455025688*10^-19 u2^3 +
4.4832557169*10^-20 u2^4) *)

nsols2 = NSolve[det2 == 0, u2]

(* Out[81]= {{u2 -> 2.9639200001}, {u2 -> 2.9639200001}, {u2 ->
0.944568025327}, {u2 -> 0.944568025327}} *)

det2 /. nsols2

(* Out[82]= {9.76746991867*10^-31, 9.76746991867*10^-31,
2.27128958953*10^-29, 2.27128958953*10^-29} *)


Is the result valid? Yes. Is it in some way better than the last one? No. The residual test is simply not showing that, despite the fact that the residuals are four orders of magnitude smaller. What is going wrong with the reasoning? It's the fact that the determinant coefficients are scaled to the square of the (already small) coefficients of the original polynomial. This again gives an indication of why one should scale input before checking a residual (scaling of the result is a somewhat more subtle matter, but also important).

If I use Solve instead of NSolve for the Kronecker product matrix determinant then I do get those smallish imaginary parts. And residuals remain small. I would not make much of that because, as noted, this residual computation is not reliably testing accuracy of solutions.

As for "exactness" of the results, one comment points out that they are in no sense exact. They are decimal approximations, suitable for an approximate input. One might or might not be able to get "better" solutions by artificially raising precision of the inputs. NSolve will do things like that under the hood, in order to either get a viable result or at least correctly assess a result for correctness via (scaled) residual computations.

Here are some other things to consider in a computation like this. One is that it might be better to replace u2 by its reciprocal so that there are no denominators to clear, then solve for the reciprocal, and invert it. The other related but more general consideration is that any time one has a mix of approximate coefficients and rational function algebra, there is room for numerical trouble. The less massaging needed, the better. In this case that amounts, again, to a problem specification that might avoid the need to clear denominators.

-

@Kellen Myers comment is useful. Since your coefficients are reals you have an approximate solution

Whenever a number carries a decimal point as your solution output, it is an approximate real.

0.//Head

(* out *)
Real


As mentioned in documentation center for Real you may change an approximate real number in an exact rational number by

approximateRealMatrix =
smatrix = {{1 - 2.96392/u2,
0.0000196744/u2}, {1. - 2.96392/u2, -2.11737*10^-10 +
0.0000196746/u2}}

exactMatrix = approximateRealMatrix /. x_Real :> Rationalize[x, 0]

sols = Det@exactMatrix // Solve[# == 0, u2] &
Det@exactMatrix /. sols[[1]]

matrix2 = KroneckerProduct[({{1, 0}, {0, 1}}), exactMatrix];
otherSol = Det@matrix2 // Solve[# == 0, u2] &
Det@exactMatrix /. otherSol[[1]]


Now you get exact results. No more decimal points indicated.

Further details on the issue may be found in the documentation center tutorial/NumericalPrecision

You mentioned, that playing with WorkingPrecision didn't change the result. Solve uses by default a WorkingPrecision value of Infinity.

Options[Solve, WorkingPrecision]

(* out *)
{WorkingPrecision -> \[Infinity]}

-
Do you really think Rationalize will turn any floating real number into some exact number -- in any context? – Kellen Myers Aug 26 '14 at 20:30
@Kellen Myers Yes, Rationalize can be used to change approximate to exact numbers. Whether this is a good thing to do, or how to do it for that matter, will depend on the underlying problem at hand. – Daniel Lichtblau Aug 26 '14 at 20:44
I believe the key word is context. Considering we start with floating point numbers, simply Rationalize-ing them will not magically change the nature of those figures or make the answer any more correct or valid, it will simply allow us to treat them as if they are. Considering hieron is only working with the six digits of precision that is displayed by default (as provided by the question), rather than the entire floating point number... what would you expect? – Kellen Myers Aug 27 '14 at 2:19