How to fit a function to 2D array when some parts of the underlaying function are known?

I have an unknown function Z=f(X,Y) where a subset of Z-values are known for a certain range of X and Y values.

Consider the following example

clc; clear; close all;
Z = [6.63900000000000e-09	3.92700000000000e-09	3.70100000000000e-09	3.37800000000000e-09	2.92600000000000e-09	2.57500000000000e-09	2.29900000000000e-09	1.89100000000000e-09	1.60700000000000e-09	1.39700000000000e-09	1.23600000000000e-09
1.13700000000000e-08	6.48300000000000e-09	6.10300000000000e-09	5.67100000000000e-09	5.02700000000000e-09	4.50200000000000e-09	4.07400000000000e-09	3.42200000000000e-09	2.95000000000000e-09	2.59200000000000e-09	2.31200000000000e-09
1.53100000000000e-08	8.55100000000000e-09	8.04300000000000e-09	7.54900000000000e-09	6.78500000000000e-09	6.14200000000000e-09	5.60700000000000e-09	4.77300000000000e-09	4.15500000000000e-09	3.67900000000000e-09	3.30000000000000e-09
1.87700000000000e-08	1.03400000000000e-08	9.72000000000000e-09	9.18400000000000e-09	8.33400000000000e-09	7.60300000000000e-09	6.98600000000000e-09	6.00700000000000e-09	5.26700000000000e-09	4.68900000000000e-09	4.22700000000000e-09
3.42700000000000e-08	1.83900000000000e-08	1.72600000000000e-08	1.66000000000000e-08	1.54900000000000e-08	1.44700000000000e-08	1.35600000000000e-08	1.20400000000000e-08	1.08200000000000e-08	9.83100000000000e-09	9.00500000000000e-09
5.17500000000000e-08	2.77800000000000e-08	2.60800000000000e-08	2.53000000000000e-08	2.40100000000000e-08	2.27700000000000e-08	2.16400000000000e-08	1.96600000000000e-08	1.80100000000000e-08	1.66200000000000e-08	1.54200000000000e-08
7.59600000000000e-08	4.14600000000000e-08	3.89400000000000e-08	3.79300000000000e-08	3.64600000000000e-08	3.50100000000000e-08	3.36600000000000e-08	3.12100000000000e-08	2.91000000000000e-08	2.72400000000000e-08	2.56100000000000e-08
1.23700000000000e-07	6.97800000000000e-08	6.54700000000000e-08	6.35700000000000e-08	6.17000000000000e-08	5.99600000000000e-08	5.83000000000000e-08	5.52400000000000e-08	5.24900000000000e-08	4.99900000000000e-08	4.77300000000000e-08
1.78100000000000e-07	1.03100000000000e-07	9.65000000000000e-08	9.27700000000000e-08	9.02000000000000e-08	8.80900000000000e-08	8.61500000000000e-08	8.25800000000000e-08	7.93200000000000e-08	7.63200000000000e-08	7.35500000000000e-08];

X = 1e-6*[0.5 0.75 1 2.5 5 7.5 10 15 20 25 30];
Y = [2.50000000000000e-07 5.00000000000000e-07 7.50000000000000e-07 1.00000000000000e-06 2.50000000000000e-06 5.00000000000000e-06 1.00000000000000e-05 2.50000000000000e-05 5.00000000000000e-05];
surf(X,Y,Z)

I want a function that predicts Z-values within the given range of X and Y values, and I know the following structure of the underlying function:

• Z is proportional to where c is an unknown constant.
• Z is proportional to where and are unknown constants.

How can I find the function f(X,Y) with Matlab?

Answer 1

[Alan Weiss]

Well, it depends what you mean by "Z is proportional to c/X and Z is proportional to c2 + c3Y + c4Y^2". You could mean that Z is proportional to c/X + c2 + c3Y + c4Y^2. In that case, you can solve for the coefficients c using backslash.

Alternatively, you could mean that Z is proportional to (c1/X + c2) * (c3 + c4Y + c5Y^2). In that case, the c vector enters nonlinearly, and you need to use lsqcurvefit. Here is a script that shows both ways. I took the liberty of scaling things to be more order one, but feel free to unscale.

Z = [6.63900000000000e-09	3.92700000000000e-09	3.70100000000000e-09	3.37800000000000e-09	2.92600000000000e-09	2.57500000000000e-09	2.29900000000000e-09	1.89100000000000e-09	1.60700000000000e-09	1.39700000000000e-09	1.23600000000000e-09
1.13700000000000e-08	6.48300000000000e-09	6.10300000000000e-09	5.67100000000000e-09	5.02700000000000e-09	4.50200000000000e-09	4.07400000000000e-09	3.42200000000000e-09	2.95000000000000e-09	2.59200000000000e-09	2.31200000000000e-09
1.53100000000000e-08	8.55100000000000e-09	8.04300000000000e-09	7.54900000000000e-09	6.78500000000000e-09	6.14200000000000e-09	5.60700000000000e-09	4.77300000000000e-09	4.15500000000000e-09	3.67900000000000e-09	3.30000000000000e-09
1.87700000000000e-08	1.03400000000000e-08	9.72000000000000e-09	9.18400000000000e-09	8.33400000000000e-09	7.60300000000000e-09	6.98600000000000e-09	6.00700000000000e-09	5.26700000000000e-09	4.68900000000000e-09	4.22700000000000e-09
3.42700000000000e-08	1.83900000000000e-08	1.72600000000000e-08	1.66000000000000e-08	1.54900000000000e-08	1.44700000000000e-08	1.35600000000000e-08	1.20400000000000e-08	1.08200000000000e-08	9.83100000000000e-09	9.00500000000000e-09
5.17500000000000e-08	2.77800000000000e-08	2.60800000000000e-08	2.53000000000000e-08	2.40100000000000e-08	2.27700000000000e-08	2.16400000000000e-08	1.96600000000000e-08	1.80100000000000e-08	1.66200000000000e-08	1.54200000000000e-08
7.59600000000000e-08	4.14600000000000e-08	3.89400000000000e-08	3.79300000000000e-08	3.64600000000000e-08	3.50100000000000e-08	3.36600000000000e-08	3.12100000000000e-08	2.91000000000000e-08	2.72400000000000e-08	2.56100000000000e-08
1.23700000000000e-07	6.97800000000000e-08	6.54700000000000e-08	6.35700000000000e-08	6.17000000000000e-08	5.99600000000000e-08	5.83000000000000e-08	5.52400000000000e-08	5.24900000000000e-08	4.99900000000000e-08	4.77300000000000e-08
1.78100000000000e-07	1.03100000000000e-07	9.65000000000000e-08	9.27700000000000e-08	9.02000000000000e-08	8.80900000000000e-08	8.61500000000000e-08	8.25800000000000e-08	7.93200000000000e-08	7.63200000000000e-08	7.35500000000000e-08];
Z = Z*1e8;
X = 1e-6*[0.5 0.75 1 2.5 5 7.5 10 15 20 25 30];
X = X*1e6;
Y = [2.50000000000000e-07 5.00000000000000e-07 7.50000000000000e-07 1.00000000000000e-06 2.50000000000000e-06 5.00000000000000e-06 1.00000000000000e-05 2.50000000000000e-05 5.00000000000000e-05];
Y = Y*1e6;
[xx,yy] = meshgrid(X,Y);
surf(xx,yy,Z)

z1 = Z(:);
x1 = xx(:);
y1 = yy(:);
v1 = [1./x1 ones(size(x1)) y1 y1.^2]; % linear problem Z ~ c1/X + c2 + c3*Y + c4*Y^2
r1 = v1\z1;

zout = v1*r1; % Implied response
zout = reshape(zout,size(Z)); % Get ready to plot
figure
surf(xx,yy,zout)

xdata = [x1 y1]; % Nonlinear response problem
rng default
x0 = randn(1,5);
options = optimoptions("lsqcurvefit","MaxFunctionEvaluations",5e5,"MaxIterations",1e5);
cout = lsqcurvefit(@twofun,x0,xdata,z1,[],[],options); % Find nonlinear parameters

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

zout2 = twofun(cout,xdata); % Get implied response
zout2 = reshape(zout2,size(Z)); % Get ready to plot
figure
surf(xx,yy,zout2)

function r = twofun(c,xdata)
x = xdata(:,1);
y = xdata(:,2);
r = (c(1)./x + c(2)) .* ( c(3) + c(4)*y + c(5)*y.^2 );
end

Alan Weiss

MATLAB mathematical toolbox documentation

[Dominik Hiltbrunner]

Hello Alan

Thank you very much for your response.

Sorry for the confusion about the "is proportional to" part. Here is the clarification:

The function f can be approximated as

Thus, your assumption for the nonlinear approximation is almost true.

I also tried to solve the problem using a least-squares approximation, but the resulting approximation is rather poor. Here is my initial attempt to solve the problem (using the normal equation):

clc; clear; close all;
Z = [6.63900000000000e-09	3.92700000000000e-09	3.70100000000000e-09	3.37800000000000e-09	2.92600000000000e-09	2.57500000000000e-09	2.29900000000000e-09	1.89100000000000e-09	1.60700000000000e-09	1.39700000000000e-09	1.23600000000000e-09
1.13700000000000e-08	6.48300000000000e-09	6.10300000000000e-09	5.67100000000000e-09	5.02700000000000e-09	4.50200000000000e-09	4.07400000000000e-09	3.42200000000000e-09	2.95000000000000e-09	2.59200000000000e-09	2.31200000000000e-09
1.53100000000000e-08	8.55100000000000e-09	8.04300000000000e-09	7.54900000000000e-09	6.78500000000000e-09	6.14200000000000e-09	5.60700000000000e-09	4.77300000000000e-09	4.15500000000000e-09	3.67900000000000e-09	3.30000000000000e-09
1.87700000000000e-08	1.03400000000000e-08	9.72000000000000e-09	9.18400000000000e-09	8.33400000000000e-09	7.60300000000000e-09	6.98600000000000e-09	6.00700000000000e-09	5.26700000000000e-09	4.68900000000000e-09	4.22700000000000e-09
3.42700000000000e-08	1.83900000000000e-08	1.72600000000000e-08	1.66000000000000e-08	1.54900000000000e-08	1.44700000000000e-08	1.35600000000000e-08	1.20400000000000e-08	1.08200000000000e-08	9.83100000000000e-09	9.00500000000000e-09
5.17500000000000e-08	2.77800000000000e-08	2.60800000000000e-08	2.53000000000000e-08	2.40100000000000e-08	2.27700000000000e-08	2.16400000000000e-08	1.96600000000000e-08	1.80100000000000e-08	1.66200000000000e-08	1.54200000000000e-08
7.59600000000000e-08	4.14600000000000e-08	3.89400000000000e-08	3.79300000000000e-08	3.64600000000000e-08	3.50100000000000e-08	3.36600000000000e-08	3.12100000000000e-08	2.91000000000000e-08	2.72400000000000e-08	2.56100000000000e-08
1.23700000000000e-07	6.97800000000000e-08	6.54700000000000e-08	6.35700000000000e-08	6.17000000000000e-08	5.99600000000000e-08	5.83000000000000e-08	5.52400000000000e-08	5.24900000000000e-08	4.99900000000000e-08	4.77300000000000e-08
1.78100000000000e-07	1.03100000000000e-07	9.65000000000000e-08	9.27700000000000e-08	9.02000000000000e-08	8.80900000000000e-08	8.61500000000000e-08	8.25800000000000e-08	7.93200000000000e-08	7.63200000000000e-08	7.35500000000000e-08];
x_range = 1e-6*[0.5 0.75 1 2.5 5 7.5 10 15 20 25 30];
y_range = [2.50000000000000e-07 5.00000000000000e-07 7.50000000000000e-07 1.00000000000000e-06 2.50000000000000e-06 5.00000000000000e-06 1.00000000000000e-05 2.50000000000000e-05 5.00000000000000e-05];
[X,Y] = meshgrid(x_range,y_range);

% Least Squares Approximation

% Convert to 1D arrays
zz = reshape(Z,1,[]);
xx = reshape(X,1,[]);
yy = reshape(Y,1,[]);

% Basis functions: {1/x, y/x, y^2/x}
syms x y
g = [1/x,y/x,y^2/x];

N = length(zz)-1;
m = length(g)-1;

% Construct design matrix
G = zeros(N+1,m+1);
for j=1:(N+1)
    for i=1:(m+1)
        G(j,i) = double(subs(g(i),{x,y},[xx(j),yy(j)]));
    end
end

% Solve normal equation
a = inv(G'*G)*G'*zz'

% Solve normal equation with QR-matrix decomposition
[Q,R] = qr(G);
Qy = Q'*zz';
a = inv(R(1:m+1,:))*Qy(1:m+1,:)

% Solve normal equation with singular-value decomposition
[U,D,V] = svd(G,'econ');
a = inv(D*V')*U'*zz'

% Verification
% Z(x,y) = a(1)*g(1) + a(2)*g(2) + ... + a(end)*g(end)
Z_approx = a(1)./X + a(2).*Y./X + a(3).*(Y.^2)./X;
surf(X,Y,Z,'FaceAlpha',0.5)     % The true values are slightly transparent
hold on
surf(X,Y,Z_approx)
e_rel_percent = 100 * abs((Z-Z_approx)./Z)

As one can see, the relative errors are quite large.

I didn't know about the lsqcurvefit function yet. According to the documentation, this function also solves the least-squares problem, so the solutions should be similar. I adjusted your solution to:

x0 = randn(1,4);

and

function r = twofun(c,xdata)
    x = xdata(:,1);
    y = xdata(:,2);
    r = (c(1)./x) .* ( c(2) + c(3)*y + c(4)*y.^2 );
end

But the resulting approximation is still rather poor. Any idea why?

[Alan Weiss]

Sorry about the poor fit. I think that the proposed model is simply not very close to the data for any values of the parameters. For example, looking at the 1-D x data for the largest value of y does not make a curve that looks like c/x to my eye. You might be able to come up with a model function that fits the data better, but you probably have good theoretical reasons for your proposed model.

Alan Weiss

MATLAB mathematical toolbox documentation

[Dominik Hiltbrunner]

There is indeed a strong theoretical reason why the Z-values are proportional to 1/x, and a not-so-strong reason for the quadratic dependency on y.

But I conclude here that the provided code works correctly; it are the model functions that need re-work.

Thank you again for your effort.

[Alex Sha]

Hi, Dominik, the fitting function below should be much better, and simple enough:

z = p1*power(y,p2+p3/x)+p4

Root of Mean Square Error (RMSE): 4.96676186258674E-9
Sum of Squared Residual: 2.44220361656497E-15
Correlation Coef. (R): 0.988334476564506
R-Square: 0.976805037566037

Parameter	Best Estimate
----------	-------------
p1	3.00671529372895E-5
p2	0.59991103826544
p3	-0.0367412799143221
p4	-1.63336827935451E-9