Fragility of matlab's solve routine, possibly a bug?

I wrote a problem set for first time matlab programmers, which involved using the solve routine. Unfortunately, they discovered the following problem, which did not endear them to either to me or matlab.

In R2016a, the solve command has the newly acquired facility to generate a numeric solution when it can't find a symbolic one, which is a very nice feature, except that it appears to be extremely fragile, treating essentially identical problems quite differently. In the example below, we compute a derivative manually, and also using matlab's symbolic diff routine. The two derivatives are identical except for rounding error of order 1e-16. In one case, the solve routine computes the correct numeric solution. In the other, it generates 17 symbolic solutions, which appear to us to be nonsense. Is this intentional behavior? Or a bug?

%%set parameters
xunder = 3;
T = 15;
b = 0.9;
r = 0.2;
rho = 2.125;
kappa = 2;
syms x nu t;
%%create base functions
c_t             = @(nu,t) nu * ((1+r) * b)^t * (1+r) * (1-b);
x_t             = @(nu,t) (1+r)^t * (xunder - nu*(1-b^t));
price           = @(x) (x/kappa)^(-1/rho);
scrapVal        = @(x) price(x) * x;
%%differentiate the scrapValue function
%  using matlab
scrapPrimeMatlab  = matlabFunction(diff(scrapVal(x)),'Vars',x);
% by hand
scrapPrimeManual  = @(x) (1 - 1/rho)*(kappa./x).^(1/rho);
% show they are essentially equal
scrapPrimeDiff  = @(x) scrapPrimeMatlab(x) - scrapPrimeManual(x);
X               = linspace(2,10,100);
disp(['Maximum difference between matlab''s deriv and manual deriv is ' num2str(max(abs(scrapPrimeDiff(X)))) ]);
disp(' ');
%%Create Control/State Functions
eqControlManual = @(x,nu) b * scrapPrimeManual(x) -  1/(c_t(nu,T));
eqControlMatlab = @(x,nu) b * scrapPrimeMatlab(x) -  1/(c_t(nu,T));
eqState         = @(x,nu) x_t(nu,T+1)-x;
% output solutions
solnManual      = solve([eqControlManual(x,nu)==0, eqState(x,nu)==0]);
disp(['Using the manual derivative, the solve command generates the numeric solution:  x=' num2str(eval(solnManual.x)) '; nu=' num2str(eval(solnManual.nu)) ]);
disp(' '); 
disp('Now running the solve command, replacing the manual derivative with matlab''s derivative')
solnMatlab      = solve([eqControlMatlab(x,nu)==0, eqState(x,nu)==0]);
solnMatlab
% different solutions, same objective function only differ by taking 
% derivative by hand with scrapPrimeManual and with matlabFunction(diff(.)) 
% where the derivative by hand and matlab differentiation are essentially
% equal

Answer 1

[John D'Errico]

Not a bug. Part of your problem is in misunderstanding floating point arithmetic, where the number 0.9 is NOT represented exactly. Nor, for example, is the number 2.5e1./1.7e1.

So lets look to see where the difference arises.

scrapPrimeMatlab
scrapPrimeMatlab = 
    @(x)x.*1.0./(x.*(1.0./2.0)).^(2.5e1./1.7e1).*(-4.0./1.7e1)+1.0./(x.*(1.0./2.0)).^(8.0./1.7e1)

scrapPrimeManual
scrapPrimeManual = 
    @(x)(1-1/rho)*(kappa./x).^(1/rho)

These are very different expressions. Subtle differences in those exponents down in the 16th decimal digit will be huge in terms of being able to solve the equations.

The fact is, symbolic algebra is not always intelligent algebra. Very often there will be terms that should combine or subtract out, but subtle differences down in those least significant bits will cause problems.

You would probably have been better off by defining numbers like 0.9 as a symbolic 0.9. Then MATLAB will be smart enough to know the number is EXACTLY 9/10. Too many of your computations were done in floating point arithmetic, and only at the very end do you tell MATLAB that this is intended to be a symbolic problem.

For example, suppose I tried this:

sym(0.9 - 0.3 - 0.6)
ans =
1/9007199254740992

Should you be surprised? Not really, IF you understand how a tool like MATLAB must do that computation. First it computes the result

0.9 - 0.3 - 0.6
ans =
   1.1102e-16

as a FLOATING point number. Then it passes that into sym. Oh, by the way, I wanted you to make that result symbolic. MATLAB cannot read your mind, at least not yet.

Sorry, but virtually computer code can be screwed up by someone who is sufficiently diligent at the task.

[Derek Wolfson]

Hey John,

Thanks for the response. I was one of the said students who hiccuped on this problem.

I understand your comment at a meta-level -- that small differences in functions can lead to huge differences when solving for zeros. However, I don't quite understand how one might still use Matlab's differentiation function within a program like this without running into problems.

Can you offer some insight? This might be rudimentary for you -- but I'm stunned by the behavior of Matlab in this particular instance.

[Leo Simon]

John, based on Derek's response, I decided that I accepted your answer too easily. Here's the same problem, with no floating point variables in sight. matlabFunction however, converts the rational numbers I've fed it (e.g., 0.9 is now 9/10, which I would hope matlab would treat as exactly that) back into floating point numbers, so that scrapPrimeMatlab looks just like it did before. So I did it a different way, and made all the relevant numbers rational floating point symbolic numbers. Then computed scrapPrimeAlt, same thing. So it seems as though Derek is right, that there's no way to bully matlab to take the derivative in a way that would avoid the problem we have later on.

T = 15;
b = 9/10;
r = 1/5;
rho = 2+1/8;
kappa = 2;
syms x nu t;
%%create base functions
c_t             = @(nu,t) nu * ((1+r) * b)^t * (1+r) * (1-b);
x_t             = @(nu,t) (1+r)^t * (xunder - nu*(1-b^t));
price           = @(x) (x/kappa)^(-1/rho);
scrapVal        = @(x) price(x) * x;
%%differentiate the scrapValue function
%  using matlab
scrapPrimeMatlab  = matlabFunction(diff(scrapVal(x)),'Vars',x);
scrapPrimeMatlab

b                   = sym(9/10);
r                   = sym(1/5);
rho                 = sym(2 + 1/8);
kappa               = sym(2);
price               = @(x) (x/kappa)^(-1/rho);
scrapVal            = @(x) price(x) * x;
scrapPrimeAlt       = matlabFunction(diff(scrapVal(x)),'Vars',x);
scrapPrimeAlt

[John D'Errico]

Hi Derek,

It is good that you are asking questions, wanting to understand.

I think your first problem MAY be in the use of functions too much here. Expressions that employ symbolic variables are already functions. For example, what if we just use syms all the way? Integers we can leave alone.

xunder = 3;
T = 15;
b = sym('0.9');
r = sym('0.2');

Note that 2.125 is indeed exactly representable as a double, so I could probably get away with leaving r as a double. But it is better to be safe than sorry.

rho = sym('2.125');
kappa = 2;
syms x nu t;
%%create base functions
c_t             = nu * ((1+r) * b)^t * (1+r) * (1-b);
x_t             = (1+r)^t * (xunder - nu*(1-b^t));
price           = (x/kappa)^(-1/rho);
scrapVal        = price * x;

Lets see where that got us:

scrapVal
scrapVal =
x/(x/2)^0.47058823529411764705882352941176

Now, lets differentiate. (Computers can be dumb, so don't be surprised at this next pair of results.)

scrapPrime = diff(simplify(scrapVal))
scrapPrime =
0.73359229710742684336247382121409/x^0.47058823529411764705882352941176

scrapPrime2 = simplify(diff(scrapVal))
scrapPrime2 =
1.3856743389806951485735616622933/x^0.47058823529411764705882352941176 - 0.65208204187326830521108784107919/x^0.47058823529411764705882352941176

Yes, you and I can see the two results are the same. And with some effort, I can probably get MATLAB to combine the latter form. Not worth the brain sweat though. :)

eqControl = b*scrapPrime - 1/subs(c_t,t,T)
eqControl =
0.66023306739668415902622643909268/x^0.47058823529411764705882352941176 - 2.6270142080490851321549115216431/nu

eqState   = subs(x_t,t,T+1)-x
eqState =
55.4652776685109248 - x - 15.062483246169511161754270763581*nu

Now just use solve.

soln      = solve([eqControl==0, eqState==0])
soln = 
    nu: [1x1 sym]
     x: [1x1 sym]

soln.nu
ans =
3.6277902362881364768614086810024

soln.x
ans =
0.82174801380353695115956262208444

In fact, I'd expect there to be no analytical solution. Pseudo-polynomial problems with fractional exponents like this almost never have an analytical solution.

Long ago I learned to go slow. Don't just throw a large block of code together and expect it to work and to do exactly as you wish, without thought on your part. Check EVERY line. Here, I don't just mean that it is mathematically correct. Does the result make sense? If symbolic,is it as simple as possible? (At least do so until you have the experience to know what you can get away with.) The point is, check each result for sanity.(I forgot to do that far too often in my early years. Ok, one truly notable time was enough to teach me that lesson.) Look to see if you could do something to improve it. For example, remember the subtle difference between simplify(diff)) and diff(simplify()). Computers can be just plain dumb sometimes. Yeah, some of the time, it is me being dumb. But not always.

As I wrote it above, see that I never had to resort to function handles. Those expressions are already functions of symbolic variables. But as I wrote them, see that by living as much as possible in the symbolic domain, MATLAB was happier. A happy computer is a good computer.