Least squares linear regression with constraints

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!

Answer 1

[Torsten]

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

[Ace]

Thank you I'll go through your solution.

[Alex Sha]

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   

Answer 2

[Bruno Luong]

"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.

[Ace]

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]

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.