Solving complicated system of symbolic equations

Hello,

I am attempting to solve a system of 3 equations that are all pretty complex and dont simplify down neatly to be used with linsolve.

The three equations at hand are:

% Given
R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];

% Create system of equations 
syms x1 x2 x3 u v 
x_vector = [x1; x2; x3];
eqn1 = x_vector'*R*x_vector == 1;
eqn2 = -(Q+Q'+2*u*R)^(-1)*(c+v*e) == x_vector;
eqn3 = (-1-e'*((Q+Q'+2*u*R)^-1)*c)/(e'*((Q+Q'+2*u*R)^-1)*e) == v;

We are given R, Q, c and e, but need to find u,v and x where x has components [x1; x2; x3]. I am wondering if this is possible on matlab and if there is another solver I should attempt to use.

If the equations above are confusing, Ive added photos of them typed out normally.

Edit: I wanted to add that I tried to use solve() with my three equations as the three arguments but I get an error that the second argument must be a vector of symbolic variables

[Walter Roberson]

please confirm that those ^-1 signify matrix inverse.

[Braden Kerr]

Yes they do. I also attempted to solve it out and break it into x, u and v matrices but the math gets very complicated, so Im wondering if its even possible to solve this on matlab without further simpplification.

[Walter Roberson]

x1 = 0.78876202436533185205858121450786
x2 = 0.64934336653846101031338613748137
x3 = -0.43810539090379286237196735198923
u = -2.6423213609632692038462718068634
v = 1.0486982172346855295291151037832

and possibly another solution that I have not checked yet.

The task gets much more difficult if you do not restrict to real values

[J. Alex Lee]

Any reason to expect analytical solutions with symbolics?

Here are 3 more real solutions (obtained numerically), according to Walter's answer it is the complete set

x1 = -0.177337903314
x2 = 0.370352455572
x3 = 0.806985447758
 u = -1.889019590663
 v = 1.993495950432

x1 = 0.687844023858
x2 = -0.389942909133
x3 = 0.702098885193
 u = 0.244348597712
 v = 2.138523547078
 
 x1 = -0.076001268894
x2 = 0.947579094374
x3 = 0.128422174517
 u = -0.806953366403
 v = 0.695475571311
 

[Braden Kerr]

To answer your question, it was for a KKT conditions problem and our professor said we could, but it gets messy. I tried, it got messy and so in talking to him he recommended using matlab instead but I havent used matlab for such a system. Im comfortable using it for other nonlinear systems or linear systems but wasnt sure how to code this one.

Thank you, this code helped me figure it out.

Answer 1

[John D'Errico]

Let me see....

eqn1 is a scalar, so it adds one piece of information. The same is true of eqn3. eqn2 is three pieces of information. And you have 5 unknowns. So in theory, it seems there are 5 pieces of information, and 5 unknowns. But that is not the complete story.

I would note that (Q + Q' + 2*u*R) will probably be singular for 3 distinct values of u. And that alone will be an issue.

R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];

syms u v
CP = det(Q + Q' + 2*u*R)

u_sing = solve(CP)
u_sing =
                         -1
- (10*181^(1/2))/59 - 20/59
  (10*181^(1/2))/59 - 20/59

So you will need to avoid those places. In fact, I would choose to restrict the search to the 4 intervals defined by the breakpoints:

uintervals = [-inf;sort(double(u_sing));inf]

uintervals =
                      -Inf
         -2.61925831306334
                        -1
          1.94129221136843
                       Inf

So how would I try to solve this?

Pick ANY value of u in one of the intervals. For example, u = 1. u lies in the interval [-1, 1.9413].

As a function of u, we can solve for v. Everything else is absolutely known.

u = 1;
% As long as we avoid the bad places for u in u_sing, M is well-posed.
M = inv(Q + Q' + 2*u*R);
v = -(1 + e'*M*c)/(e'*M*e);
x = -M*(c + v*e);

% Does x satisfy the requirement that 
%  x'*R*x == 1?
x'*R*x 
         0.276636817658454

So, no. But we can view this as a problem of a ONE variable search, where we change u in each of the intervals we have identified. Find the value of u such that x'*R*x == 1, or is as close as possible to that goal.

function eqn1resid = uobj(u,Q,R,e,c)
  M = inv(Q + Q' + 2*u*R);
  v = -(1 + e'*M*c)/(e'*M*e);
  x = -M*(c + v*e);
  eqn1resid = abs(x'*R*x - 1);
end

So we wish to MINIMIZE eqn1resid, as a function of u.

[uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),-10,uintervals(2))
uopt =
         -2.64230839840025
fval =
      5.88144443320893e-05
exitflag =
     1
     
     [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(2),uintervals(3))
uopt =
         -1.88901086809246
fval =
       3.4982453337884e-05
exitflag =
     1
     
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),uintervals(4))
uopt =
         0.244333979030155
fval =
      4.84817244350566e-05
exitflag =
     1

Finally, in the topmost interval for u, fminbnd does not find any solution.

>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(4),100)
uopt =
          1.94134412859712
fval =
         0.823090976327026
exitflag =
     1

The point being in all of this, is you do not need to do anything fancy. Double precision arithmetic seems to work.

Will you find an algebraic solution? Probably not. But you should be able to find numerical solutions.

[J. Alex Lee]

neat, didn't quite see that it can be posed as 1D in one unknown u

But, why do you want to minimize abs(x'Rx-1) versus zero x'Rx-1 with fzero (or a bracketing algorithm)?

[John D'Errico]

An excellent question, one that I considered myself in making the choice of solver.

Note that it appears there is no solution in the topmost interval for u. See that fminbnd returns a number that is not effectively zero. I'd probably want to do some more investigation to convince myself there is truly no solution up there. But fminbnd seems to think there is no solution. fzero might come to the same conclusion, but there is a more serious problem with fzero.

Fzero can operate without a bracket, if I give it just a single start point. But if I do not provide a bracket for fzero, I cannot insure it will find a solution in the interval I want to search. And if I do provide a bracket for fzero, then fzero will fail if the objective has the same sign on each end of the bracket. In any case, since I do wish to solve the problem on a bracketed interval, then fminbnd seems the perfectly natural solver.

Essentially, fminbnd becomes a somewhat more robust solution. You just need to verify that fval is approximately zero at the solution point.

[J. Alex Lee]

Good points. I guess still, depending on the minimization algorithm, you wouldn't necessarily be able to find all roots within a particular interval if there are multiple ones. I guess case in point, this method missed the root u = -0.80695 in the 3rd interval.

I'm fascinated that the system is reducible to a polynomial in u, as Walter showed with Maple, which must be the better answer for this specific problem...but numerically (if we didn't know the equation for u was polynomial with determinable coefficients) I guess one can scan over some reasonable range of u to look for cross-overs, especially if some more analysis of the equations and known values could narrow the range from [-Inf,Inf] to something finite.

Any way to do better than the [-Inf,Inf] outer bracket?

[John D'Errico]

That Walter was able to solve the problem using symbolic methods is a testament to the power of those methods. It also points out there are almost always multiple ways to solve any given problem.

I was lazy yesterday when I wrote my answer though. One of the things I should probably have done was to plot the resulting function I was optimizing, thus plotting x'*R*x-1, as a function of u. That is my fault, because it fails one of the most important rules I tell people to do: PLOT EVERYTHING. And then plot something else. Plot it differently. Think about what you see. That plot forces the analyst to verify they have done something reasonable. Does what you see make sense?

The other thing a plot does is to help you to visualize what is happening. Is there a root on some interval? Are there multiple roots? Is your problem practically solvable?

I'll concede that one thing such an intuitive graphical approach might do is to target the user to the use of numerical solvers, but even with symbolic methods, you should still follow the rule that if you can plot what you are doing, then you should do so, before you try to solve anything.

function res = xRxm1(u,Q,R,e,c)
  res = zeros(size(u));
  for i = 1:numel(u)
    M = inv(Q + Q' + 2*u(i)*R);
    v = -(1 + e'*M*c)/(e'*M*e);
    x = -M*(c + v*e);
    res(i) = x'*R*x - 1;
  end
end

I've removed the abs there, so we can look at what it produces more intelligently in terms of finding a zero. I also vectorized the code, so it can handle multiple values for u at once.

R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
c = [-1; -1; -1];
e = [1; 1; 1];

As I said, the intervals of interest for u should be defined by the breaks:

uintervals = [-Inf, -2.61925831306334, -1, 1.94129221136843, Inf];

The inner breaks occur where our matrix was singular, as a function of u. So this method must fail at those values for u. We can verify that fact.

>> xRxm1(-1,Q,R,e,c)
Warning: Matrix is singular to working precision. 
> In xRxm1 (line 4)
 
ans =
   NaN

Oh - there is no outer interval to consider. The first interval goes from [-inf,-2.6]. I reduced +/-inf to +/-10 in the code I played with before, because numerical solvers tend not to like infs. So the top interval is [1.94,inf].

Now, plot this function over an interval of interest. When we do so, we will see something interesting. Remember, in the plot here, we are looking for a zero crossing.

fplot(@(u) xRxm1(u,Q,R,e,c),[-.9,1.94])

So it looks like the function has a singularity around x=-0.28 or so. In fact, it almost looks like no zero crossing exists. But if we break the problem down into two further sub-intervals, we see this:

fplot(@(u) xRxm1(u,Q,R,e,c),[-.1,1.94])
yline(0);

fplot(@(u) xRxm1(u,Q,R,e,c),[-.999,-0.4])
yline(0);

So as you see, fzero would have failed to produce a root, if we just gave it the complete interval [-1,1.94]. That is because an objective function that fzero would use will produce end points that have the same sign at each end of that interval. Whereas fminbnd does produce a solution. It will find ONE of those roots.

>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),uintervals(4))
uopt =
         0.244333979030155
fval =
      4.84817244350566e-05
exitflag =
     1

And now I'll force it to find the other:

     
>> [uopt,fval,exitflag] = fminbnd(@(u) uobj(u,Q,R,e,c),uintervals(3),-.4)
uopt =
        -0.806937125898144
fval =
      4.91463840461837e-05
exitflag =
     1

I could still use fzero. But I need to be careful about the bracket choice.

[uobj,fval,exitflag] = fzero(@(u) xRxm1(u,Q,R,e,c),[-.999,-.4])
uobj =
        -0.806953366404047
fval =
     -2.22044604925031e-16
exitflag =
     1

Anyway, if we look at each of those intervals, I'd bet that there may actually be TWO such real solutions in most of the intervals I have identified. Except for the top most interval, where none seems to exist at all.

The only other point of interest is to identify why there is a singularity in the function I have been plotting. But really, I need to stop somewhere. :)

As I said at the beginning of thos coment though, I really should have plotted the thing first. That would have made my answer more complete. And vastly longer. :(

Answer 2

[Walter Roberson]

Maple gives a fast response that

 RR = roots([3353, 17080, 24623, 6420, - 3300])
 u = RR
 v = 10059/11420*RR^3 + 3843/1142*RR^2 + 32129/11420*RR + 1413/1142
 x1 = -466067/2969200*RR^3 - 26173/296920*RR^2 + 2249503/2969200*RR + 151511/296920
 x2 = -55085/118768*RR^3 - 117787/59384*RR^2 - 1787385/831376*RR + 108319/415688
 x3 = 17723/28550*RR^3 + 11829/5710*RR^2 + 278251/199850*RR + 9159/39970

There are analytic solutions for RR, as it is a degree 4 polynomial. The numeric values are approximately [0.2443485977827867, -0.8069533664040474, -1.889019590665992, -2.642321360963269] . The analytic solutions are not worth posting as they are too long.

[J. Alex Lee]

so, maple just spat out that the problem reduces a 4th order polynomial in u?

and also, that catches a mistake of one of my invalid solutions since you show there should only be 4 real ones

[Walter Roberson]

I was surprised how quickly Maple solved it; I had already spent about half an hour trying to get a solution out of MATLAB by analyzing the circumstances under which the ReturnConditions could be true, and Maple responded with a complete solution almost immediately.

Maple does not return in exactly the form I posted: each variable is in terms of RootOf() expressions, occuring multiple times for some of the elements. Pulling those RootOf out into a parameter is easy though:

SOL := solve(EQN);
R := indets(SOL, RootOf)[1]; 
S := subs(R=RR, SOL);

Answer 3

[J. Alex Lee]

I like John's comment about visualizing. To expand on this, and Walter's comment that it is trickier to do this in the complex plane, here is a numerical exploration for complex solutions assuming is the conjugate transpose

clc;
close all;
clear

R = [1 -.5 -.4; -.5 1 0; -.4 0 1];
Q = [1 0 -2; 0 1 0; -2 0 1];
cnt = [-1; -1; -1];
e = [1; 1; 1];

% real roots by Walter (for visualization)
uRoots = roots([3353, 17080, 24623, 6420, - 3300])

% singular points in u by John (for visualization
uSing = [-1, - (10*181^(1/2))/59 - 20/59, (10*181^(1/2))/59 - 20/59]

% define a grid on the complex plane
[uR,uI] = ndgrid(linspace(-3,2,200),linspace(-0.8,0.8,200));
resMat = nan(size(uR));

% calculate residual for each grid point
for k = 1:numel(uR)
	u = uR(k) + 1i*uI(k);
	[r,v,x] = resfun1(u , R,Q,cnt,e);
	% the residual should be real because of the way it is defined (?)
	resMat(k) = real(r);
end

% visualize
figure(1); cla; hold on;
% the actual contour map of residual value
contourf(uR,uI,resMat,linspace(-1,1,35)*0.8,'LineStyle','none')
% draw the curves of solution
contour(uR,uI,resMat,[0,0],'EdgeColor','w','LineWidth',1)
colormap turbo
shading interp
colorbar
set(gca,'DataAspectRatio',[1 1 1])

% plot the singular points on real line
plot(uSing,zeros(size(uSing)),'xw','LineWidth',1,'MarkerSize',6)

% plot the real roots on the real line
plot(uRoots,zeros(size(uRoots)),'ow','LineWidth',1,'MarkerSize',4)

% plot the centers of ostensible circles
cnt = [mean(uRoots(1:2)),mean(uRoots(3:4))]
rad = (uRoots([2,4])'-c)
plot(cnt,zeros(size(cnt)),'+k')

xlabel('real part of $u$','Interpreter','latex')
ylabel('imaginary part of $u$','Interpreter','latex')
title('residual $\left(1-x^T R x\right)$','Interpreter','latex')

function [r,v,x] = resfun1(u , R,Q,c,e)

Z = Q+Q'+2*u*R;
Zi = Z\eye(3);

v = -(1+e'*Zi*c)/(e'*Zi*e);
x = -Zi*(c+v*e);

r = 1 - x'*R*x; % conjugate transpose ensures real result (I think)

end

The solution looks to be circular rings about the mid-points of the 2 consecutive pairs of real solutions found by Walter. The midpoints are the singularities in the residuals found by John.

So there may be a bit nicer structure in the solutions if analyzed more deeply?