How do I solve over determined system of non linear equation?
4 views (last 30 days)
Show older comments
Ahmad Adee
on 26 May 2018
Commented: John D'Errico
on 28 Feb 2021
Suppose I have three matrices A,B,X. A and B being now and size 200x4. X being unknown with size 3x4. The system is AX=B.
syms x1 x2 x3 x4 y1 y2 y3 y4 z1 z2 z3 z4
A = rand(200,4);
X = [ x1 x2 x3 x4; y1 y2 y3 y4; z1 z2 z3 z4];
B = rand(200,4);
I want solve for X while keeping x1 x2 x3 y1 y2 y3 z1 z2 z3 values constrained to [-1,1]. What is the best approach?
2 Comments
John D'Errico
on 27 May 2018
Edited: John D'Errico
on 27 May 2018
I'm not sure why you think it is nonlinear. Purely linear. Admittedly, you have a size problem, so matrix multiplication does not conform.
Do I understand that you want the last row of the matrix X to be [0 0 0 1]?
Accepted Answer
John D'Errico
on 27 May 2018
Edited: John D'Errico
on 27 May 2018
Why do people think they need to use syms, just because they don't know the value of a result? (I think that somehow comes from the design of the symbolic toolbox and how it interacts with MATLAB.) Basic linear algebra suffices here. And you CAN solve the problem. You just need to accept that since it is over-determined, there is no exact solution in general.
A = rand(200,4);
B = rand(200,4);
Now, you want to solve the problem A*X=B, where X has a fixed, known row. Or apparently it may even be multiple rows in your actual problem. So I'll show how to solve it for an array of size 4x4, but with some set of known rows.
fixedRowIndex = 4;
X_fixed = [0 0 0 1];
estRowIndex = setdiff(1:size(A,2),fixedRowIndex);
Now, we can imagine partitioning the matrix X into two sub-matrices, they will look sort of like this:
% X = [X_est;X_fixed]
A is a 200x4 array. But when we multiply it by X, the problem splits into two pieces.
% A*X = A(:,estRowIndex)*X_est + A(:,fixedRowIndex)*X_fixed
So, does this make sense? If so, then you need to recognize that the second term is fully KNOWN. So we can move it to the right hand side. The problem reduces to solving for the unknown 3x4 array Xest:
% A(:,estRowIndex)*X_est = B - A(:,fixedRowIndex)*X_fixed
Now that we have eliminated those known rows from X, solve it simply, using backslash.
X_est = A(:,estRowIndex)\(B - A(:,fixedRowIndex)*X_fixed)
X_est =
0.45678 0.3027 0.2734 0.027561
0.25445 0.39719 0.3131 0.0085818
0.23474 0.28378 0.32242 -0.067524
X(estRowIndex,:) = X_est;
X(fixedRowIndex,:) = X_fixed
X =
0.45678 0.3027 0.2734 0.027561
0.25445 0.39719 0.3131 0.0085818
0.23474 0.28378 0.32242 -0.067524
0 0 0 1
The result is the least squares solution, consistent with the constraint that X has some fixed, known in advance rows. The nice thing is as I wrote it, if you have multiple fixed, known rows, it is trivially extensible.
Could I have done this using lsqlin? Yes, of course. But why bother?
2 Comments
Ameer Hamza
on 27 May 2018
Edited: Ameer Hamza
on 27 May 2018
This solution should be fine for a non-constrained case. But the question required that element of X are constrained to [-1, 1]. For example
A = rand(200,4);
B = 10*rand(200,4);
fixedRowIndex = 4;
X_fixed = [0 0 0 1];
estRowIndex = setdiff(1:size(A,2),fixedRowIndex);
X_est = A(:,estRowIndex)\(B - A(:,fixedRowIndex)*X_fixed);
X(estRowIndex,:) = X_est;
X(fixedRowIndex,:) = X_fixed
X =
3.7302 2.2907 3.1585 3.1444
1.9967 3.4358 2.7390 3.2424
3.4536 2.6884 3.0786 1.9789
0 0 0 1.0000
The X elements no longer belong to [-1, 1].
John D'Errico
on 27 May 2018
Edited: John D'Errico
on 27 May 2018
Sigh. Still trivial. Of course, lsqlin will now be needed. Must I point out that you first wanted to "solve" the problem using X = A\B?
So how do we EFFICIENTLY solve this problem using LSQLIN and bound constraints?
LSQLIN applies to linear problems where you have bound constraints. Of course, LSQLIN will only solve one right hand side at a time. And really this problem has 4 right hand sides. backslash solved that in one line, but really, it was 4 basically parallel problems.
So the answer is one of two solutions. Simplest is to loop. This was the original solve.
X_est = A(:,estRowIndex)\(B - A(:,fixedRowIndex)*X_fixed)
But, IF we wanted to solve using LSQLIN and bound constraints, with a loop...
options = optimset('lsqlin');
options.Display = 'none';
Bhat = B - A(:,fixedRowIndex)*X_fixed;
A_est = A(:,estRowIndex);
nrhs = size(Bhat,2);
X = zeros(size(A,2),size(B,2));
X(fixedRowIndex,:) = X_fixed;
lb = repmat(-1,size(estRowIndex));
ub = repmat(1,size(estRowIndex));
for i = 1:nrhs
X(estRowIndex,i) = lsqlin(A_est,Bhat(:,i),[],[],[],[],lb,ub,[],options);
end
No use of fmincon needed.
Finally, could I have re-written this to avoid the loop completely? Well, yes, all in one call to lsqlin. The solution would be somewhat less efficient unless it is done carefully though.
More Answers (2)
Ameer Hamza
on 27 May 2018
Edited: Ameer Hamza
on 27 May 2018
Since the system is Overdetermined, the unique solution is highly unlikely. The best you can expect is the solution to the following least squares minimization problem
argmin_X mean2(abs(A*X-B))
You can get it like this
A = rand(200,4);
B = rand(200,4);
X = A\B
Edit: since you mentioned other constraints to be fulfilled, therefore you can solve the modified problem as follow
A = rand(200,4);
B = rand(200,4);
f = @(x) sum(sum((A*[x; 0 0 0 1]-B).^2)); % constraining the last row of x to be [0 0 0 1]
X = fmincon(f, zeros(3, 4), [], [], [], [], -1*ones(3, 4), ones(3, 4));
X = [X; 0 0 0 1];
This will force the last row of X to be [0 0 0 1] and also put constraints on other elements of X to be between -1 to 1.
Note the updated objective function use element-wise norm, which is same as used by MATLAB's lsqlin, thereby producing the same solution.
3 Comments
Ameer Hamza
on 27 May 2018
The objective function of this fully linear least square problem is non-linear. I agree with your comment that lsqlin() is a possible solution, but fmincon is solving the same problem in a more compact and intuitive way to solve a non-linear objective function.
John D'Errico
on 28 Feb 2021
But the point is, the OBJECTIVE function is not nonlinear. It is quadratic. Using a nonlinear solver will force the tool to differentiate the problem, and then iterate to find a solution on a problem that requires no iteration. In the end, you will get a result that will be less accurate, since this is only an iterative method that will try then to converge, but only to within a tolerance.
Walter Roberson
on 26 May 2018
That is not possible. If A is something by 4, and X is 3 by something, then A*X would be trying to do algebraic matrix multiplication between arrays of inconsistent sizes.
In order to be able to calculate A*X with A being something by 4, then X would have to be 4 by something. If the result is to be 200x4 with A being 200 x 4, then X would have to be 4 x 4, definitely not 3 x 4.
See Also
Categories
Find more on Linear Least Squares in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!