Least squares linear regression with constraints

23 views (last 30 days)
Ace
Ace on 30 Jun 2023
Moved: Bruno Luong on 4 Aug 2023
Hi everyone!
I am new to MATLAB and I have the following issue. I have the following data set:
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
where x11 and x12 are the measured inputs, and y is the known output.
I need to solve for the unknown coefficients a1, a2 and a3 in the following original equation:
y = a1*(x11^a2)*(x12^a3)
This is a non-linear problem of course, but I do linearize it by taking log on both sides, so that I have:
y2 = log10(y);
x1 = log10(x1);
x2 = log10(x2);
I then perform a linear regression analysis using the pinv function which makes use of a pseudoinverse matrix to find the coefficients.
My current code compiled based on the information I found on the Internet is following:
clear all close all
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
x1=log10(x11);
x2=log10(x12);
y2=log10(y);
n=length(x1);
a=[ones(n,1), x1', x2'];
c=pinv(a)*y2';
a1=c(1,:);
a2=c(2,:);
a3=c(3,:);
%here I take antilog to get back the original power relations:
y_pred = 10^(log10(a1)).*(x11.^a2).*(x12.^a3);
Although I do get a solution for all of the coefficients I am having two issues:
  1. I am not really satisfied with the output (poor R^2 coefficient)
  2. I want to make some constraints so that all of the coefficient stay positive, for example: using lsqlin didn't help much as the solution the is highly biased by my initial suggestions and the coefficient values either take the upper or lower bound value, which seems not to be correct.
I obtained a better solution using fminsearch approach without pre-linearizing the equation, but was curious if I can get comparable results taking the linear way.
I would be grateful for your help!

Answers (2)

Torsten
Torsten on 30 Jun 2023
Edited: Torsten on 30 Jun 2023
You made a mistake in computing the correct coefficients (see below).
clear all close all
x11 = [0.091, 0.068, 0.086, 0.091, 0.097, 0.064, 0.052, 0.066, 0.074];
x12 = [0.707, 0.920, 0.491, 0.527, 0.656, 0.474, 0.49, 0.267, 0.552];
y = [0.07, 0.029, 0.063, 0.077, 0.155, 0.043, 0.023, 0.057, 0.085];
x1=log10(x11);
x2=log10(x12);
y2=log10(y);
n=length(x1);
a=[ones(n,1), x1', x2'];
c=a\y2';
a1_lin=10^c(1,:)
a1_lin = 34.9283
a2_lin=c(2,:)
a2_lin = 2.5909
a3_lin=c(3,:)
a3_lin = -0.4904
%here I take antilog to get back the original power relations:
y_pred = a1_lin.*(x11.^a2_lin).*(x12.^a3_lin);
norm(y_pred-y)
ans = 0.0695
fun = @(x)x(1)*x11.^x(2).*x12.^x(3)-y;
sol = lsqnonlin(fun,[a1_lin a2_lin a3_lin])
Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality tolerance.
sol = 1×3
59.7121 2.7205 -0.1921
a1_nonlin = sol(1)
a1_nonlin = 59.7121
a2_nonlin = sol(2)
a2_nonlin = 2.7205
a3_nonlin = sol(3)
a3_nonlin = -0.1921
norm(fun(sol))
ans = 0.0666
  2 Comments
Ace
Ace on 30 Jun 2023
Thank you I'll go through your solution.
Alex Sha
Alex Sha on 4 Aug 2023
Moved: Bruno Luong on 4 Aug 2023
If using direct nonlinear fitting
a1 59.737732722511
a2 2.72067588034148
a3 -0.192039313150924
While doing by linear transformation:
a1 34.9283448971434
a2 2.59089956584206
a3 -0.490446537714015
There are significant difference between two methods.
If want to all positive values:
a1 52.3606801051862
a2 2.62431938697266
a3 0

Sign in to comment.


Bruno Luong
Bruno Luong on 30 Jun 2023
Edited: Bruno Luong on 30 Jun 2023
"This is a non-linear problem of course, but I do linearize it by taking log on both side"
That's the problem. When you transform problem by taking the log, if you have Gaussian noise on data, it will no longer be Gaussian, let alone normalized them with covariance matrix. So least-squares with transformed model is not correct way to sove the problem.
You need to use non-linear regression, fmincon or lsqnonlin (optimization toolbox required), that is a good way to go. fminsearch migh be not a good tool numerically, with poor convergence.
  2 Comments
Ace
Ace on 30 Jun 2023
Thanks for your suggestion. But how do I put my contraint (positive values) not to force the model to show specific outcomes? Using upper and lower bounds seems too aggressive: the result might be actually equal to one of those. Thanks.
Bruno Luong
Bruno Luong on 30 Jun 2023
If the constraints respect the physics then it is not agressive.
You have a range of ways to enter the contarints in fmincon and lsqnonlin, just read the doc.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!