Optimization method based on the Nonlinear least squares as well as fmincon for the defined objective function (Non-linear function).
Hi there,
I am currently working on a Matlab code for my task (optimization task), and I tried two different methods;
Method A: the fmincon was used to minimize the objective function in which different algorithms, such as interior-point, sqp, sqp-legacy, interior-point, and active-set, were applied.
Method B: Nonlinear least squares optimization was combined with three algorithms, levenberg-marquardt, and trust-region-reflective.
The best results I got now are (from the Nonlinear least squares optimization + ‘interior-point’):
RMSE for X, Y, Z Direction(mm) : [1.18, 0.16, 1.33] mm; however, I want to reduce it as far as it is possible.
@ A quick guide for the code, and how it works.
1. I do have 3D data sets, and I used some mathematical methods to convert these 3D data into 2D data (Step 3).
2. The goal is to use some inverse method to get the same 3D datasets; therefore, after defining the objective function, I applied different optimization methods to address this task.
Any assistance or alternative approaches would be greatly appreciated.
Thank you in advance!
% @ Payam Samadi (2025.3.15) NLSP (Nonlinear least squares optimization)
% In step 3, the 2D projection data was extracted from the 3D data at
% different projections.
% In step4, the objective function is defined.
% In step 5, the Nonlinear least squares optimization is used to minimize
% the objective function.
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Note: It is important to carefully consider the initial guesses for
% target locations, as well as the lower bounds (lb) and upper bounds (ub).
tic;
clc;
clear;
close all;
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread(‘Sample 1.xlsx’); % Data include 12 (sec)
% Load_Data = xlsread(‘Sample 2.xlsx’); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
[~,numT]=size(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha); % view angles (radians)
%% Step 3: Projection Model: Compute 2D projections
Xp=zeros(1,num_projections); % allocate array
Yp=zeros(1,num_projections);
for i=1:num_projections
f_theta = (SAD – (true_T(i,1) .* sin(theta_rad(i)) + true_T(i,3) .* cos(theta_rad(i)))) ./ SID;
Xp(1, i) = (true_T(i,1) .* cos(theta_rad(i)) – true_T(i,3) .* sin(theta_rad(i))) ./ f_theta;
Yp(1, i) = true_T(i,2) ./ f_theta;
end
xp = Xp’;
yp = Yp’;
%% Step 4: Define the objective function
% Define the objective function
objFun = @(T) computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections);
%% Step 5: Minimize the objective function using lsqnonlin
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Define bounds
lb = []; % Lower bounds
ub = []; % Upper bounds
% Define optimization options
options = optimoptions(‘lsqnonlin’, …
‘Algorithm’, ‘interior-point’, …
‘Display’, ‘iter’, …
‘MaxIterations’, 300, …
‘MaxFunctionEvaluations’, 20000, …
‘TolFun’, 1e-10, …
‘TolX’, 1e-10, …
‘StepTolerance’, 1e-8, …
‘FiniteDifferenceStepSize’, 1e-3, …
‘UseParallel’, true);
% Generate a better initial guess
T0 = ones(3, num_projections);
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
for i = 1:num_projections
objFun_i = @(T) computeProjectionError(T, xp(i), yp(i), theta_rad(i), SAD, SID, 1);
estimated_T(:, i) = lsqnonlin(objFun_i, T0(:, i), lb, ub, options);
end
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(1,:);
Estimate_3D_Y = estimated_T(2,:);
Estimate_3D_Z = estimated_T(3,:);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf([‘RMSE for X, Y, Z Direction(mm): ‘ …
‘[%.2f, %.2f, %.2f] mmn’], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc;
%% Objective function: Computes error between observed and estimated projections
function error = computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections)
for i=1:num_projections
x_est = T(1,:);
y_est = T(2,:);
z_est = T(3,:);
% Compute estimated projection using the perspective transformation
f_theta = (SAD – (x_est(i) .* sin(theta_rad(i)) + z_est(i) .* cos(theta_rad(i))))./ SID;
xp_est(i) = (x_est(i) .* cos(theta_rad(i)) – z_est(i) .* sin(theta_rad(i))) ./ f_theta;
yp_est(i) = y_est(i) ./ f_theta;
end
%% Compute residual error (Method 1)
% Compute the residual error
error_x = xp – xp_est;
error_y = yp – yp_est;
% Return residuals for least squares minimization
error = [error_x; error_y];
end
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,’omitnan’));
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure(‘WindowState’, ‘maximized’)
subplot(3,1,1);
plot(Time, Estimate_3D_X,’r’);
hold on;
plot(Time, xt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘X Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,’r’);
hold on;
plot(Time, yt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Y Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,’r’);
hold on;
plot(Time, zt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Z Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
endHi there,
I am currently working on a Matlab code for my task (optimization task), and I tried two different methods;
Method A: the fmincon was used to minimize the objective function in which different algorithms, such as interior-point, sqp, sqp-legacy, interior-point, and active-set, were applied.
Method B: Nonlinear least squares optimization was combined with three algorithms, levenberg-marquardt, and trust-region-reflective.
The best results I got now are (from the Nonlinear least squares optimization + ‘interior-point’):
RMSE for X, Y, Z Direction(mm) : [1.18, 0.16, 1.33] mm; however, I want to reduce it as far as it is possible.
@ A quick guide for the code, and how it works.
1. I do have 3D data sets, and I used some mathematical methods to convert these 3D data into 2D data (Step 3).
2. The goal is to use some inverse method to get the same 3D datasets; therefore, after defining the objective function, I applied different optimization methods to address this task.
Any assistance or alternative approaches would be greatly appreciated.
Thank you in advance!
% @ Payam Samadi (2025.3.15) NLSP (Nonlinear least squares optimization)
% In step 3, the 2D projection data was extracted from the 3D data at
% different projections.
% In step4, the objective function is defined.
% In step 5, the Nonlinear least squares optimization is used to minimize
% the objective function.
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Note: It is important to carefully consider the initial guesses for
% target locations, as well as the lower bounds (lb) and upper bounds (ub).
tic;
clc;
clear;
close all;
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread(‘Sample 1.xlsx’); % Data include 12 (sec)
% Load_Data = xlsread(‘Sample 2.xlsx’); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
[~,numT]=size(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha); % view angles (radians)
%% Step 3: Projection Model: Compute 2D projections
Xp=zeros(1,num_projections); % allocate array
Yp=zeros(1,num_projections);
for i=1:num_projections
f_theta = (SAD – (true_T(i,1) .* sin(theta_rad(i)) + true_T(i,3) .* cos(theta_rad(i)))) ./ SID;
Xp(1, i) = (true_T(i,1) .* cos(theta_rad(i)) – true_T(i,3) .* sin(theta_rad(i))) ./ f_theta;
Yp(1, i) = true_T(i,2) ./ f_theta;
end
xp = Xp’;
yp = Yp’;
%% Step 4: Define the objective function
% Define the objective function
objFun = @(T) computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections);
%% Step 5: Minimize the objective function using lsqnonlin
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Define bounds
lb = []; % Lower bounds
ub = []; % Upper bounds
% Define optimization options
options = optimoptions(‘lsqnonlin’, …
‘Algorithm’, ‘interior-point’, …
‘Display’, ‘iter’, …
‘MaxIterations’, 300, …
‘MaxFunctionEvaluations’, 20000, …
‘TolFun’, 1e-10, …
‘TolX’, 1e-10, …
‘StepTolerance’, 1e-8, …
‘FiniteDifferenceStepSize’, 1e-3, …
‘UseParallel’, true);
% Generate a better initial guess
T0 = ones(3, num_projections);
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
for i = 1:num_projections
objFun_i = @(T) computeProjectionError(T, xp(i), yp(i), theta_rad(i), SAD, SID, 1);
estimated_T(:, i) = lsqnonlin(objFun_i, T0(:, i), lb, ub, options);
end
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(1,:);
Estimate_3D_Y = estimated_T(2,:);
Estimate_3D_Z = estimated_T(3,:);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf([‘RMSE for X, Y, Z Direction(mm): ‘ …
‘[%.2f, %.2f, %.2f] mmn’], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc;
%% Objective function: Computes error between observed and estimated projections
function error = computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections)
for i=1:num_projections
x_est = T(1,:);
y_est = T(2,:);
z_est = T(3,:);
% Compute estimated projection using the perspective transformation
f_theta = (SAD – (x_est(i) .* sin(theta_rad(i)) + z_est(i) .* cos(theta_rad(i))))./ SID;
xp_est(i) = (x_est(i) .* cos(theta_rad(i)) – z_est(i) .* sin(theta_rad(i))) ./ f_theta;
yp_est(i) = y_est(i) ./ f_theta;
end
%% Compute residual error (Method 1)
% Compute the residual error
error_x = xp – xp_est;
error_y = yp – yp_est;
% Return residuals for least squares minimization
error = [error_x; error_y];
end
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,’omitnan’));
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure(‘WindowState’, ‘maximized’)
subplot(3,1,1);
plot(Time, Estimate_3D_X,’r’);
hold on;
plot(Time, xt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘X Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,’r’);
hold on;
plot(Time, yt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Y Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,’r’);
hold on;
plot(Time, zt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Z Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
end Hi there,
I am currently working on a Matlab code for my task (optimization task), and I tried two different methods;
Method A: the fmincon was used to minimize the objective function in which different algorithms, such as interior-point, sqp, sqp-legacy, interior-point, and active-set, were applied.
Method B: Nonlinear least squares optimization was combined with three algorithms, levenberg-marquardt, and trust-region-reflective.
The best results I got now are (from the Nonlinear least squares optimization + ‘interior-point’):
RMSE for X, Y, Z Direction(mm) : [1.18, 0.16, 1.33] mm; however, I want to reduce it as far as it is possible.
@ A quick guide for the code, and how it works.
1. I do have 3D data sets, and I used some mathematical methods to convert these 3D data into 2D data (Step 3).
2. The goal is to use some inverse method to get the same 3D datasets; therefore, after defining the objective function, I applied different optimization methods to address this task.
Any assistance or alternative approaches would be greatly appreciated.
Thank you in advance!
% @ Payam Samadi (2025.3.15) NLSP (Nonlinear least squares optimization)
% In step 3, the 2D projection data was extracted from the 3D data at
% different projections.
% In step4, the objective function is defined.
% In step 5, the Nonlinear least squares optimization is used to minimize
% the objective function.
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Note: It is important to carefully consider the initial guesses for
% target locations, as well as the lower bounds (lb) and upper bounds (ub).
tic;
clc;
clear;
close all;
%% Step 1: Load xls file (3D Tumor position)
Load_Data = xlsread(‘Sample 1.xlsx’); % Data include 12 (sec)
% Load_Data = xlsread(‘Sample 2.xlsx’); % Data include 60 (sec)
Time = Load_Data(:,1); % Time (s)
xt=Load_Data(:,2); % xt for targets
yt=Load_Data(:,3); % yt for targets
zt=Load_Data(:,4); % zt for targets
true_T=[xt,yt,zt]; % target locations
[~,numT]=size(true_T); % number of targets
%% Step 2: Define imaging system parameters
SAD = 100; % source-axis distance (cm)
SID = 150; % source-image plane distance (cm)
num_projections = length(Time); % number of different views
% Generate projection angles
alpha = linspace(0, 360, num_projections);
% alpha=30;
theta_rad = deg2rad(alpha); % view angles (radians)
%% Step 3: Projection Model: Compute 2D projections
Xp=zeros(1,num_projections); % allocate array
Yp=zeros(1,num_projections);
for i=1:num_projections
f_theta = (SAD – (true_T(i,1) .* sin(theta_rad(i)) + true_T(i,3) .* cos(theta_rad(i)))) ./ SID;
Xp(1, i) = (true_T(i,1) .* cos(theta_rad(i)) – true_T(i,3) .* sin(theta_rad(i))) ./ f_theta;
Yp(1, i) = true_T(i,2) ./ f_theta;
end
xp = Xp’;
yp = Yp’;
%% Step 4: Define the objective function
% Define the objective function
objFun = @(T) computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections);
%% Step 5: Minimize the objective function using lsqnonlin
% Algorithm:
% 1. ‘levenberg-marquardt’ : [1.21, 0.17, 1.39] mm,
% 2. ‘trust-region-reflective’ : [1.21, 0.17, 1.39] mm,
% 3. ‘interior-point’ : [1.18, 0.16, 1.33] mm
% Define bounds
lb = []; % Lower bounds
ub = []; % Upper bounds
% Define optimization options
options = optimoptions(‘lsqnonlin’, …
‘Algorithm’, ‘interior-point’, …
‘Display’, ‘iter’, …
‘MaxIterations’, 300, …
‘MaxFunctionEvaluations’, 20000, …
‘TolFun’, 1e-10, …
‘TolX’, 1e-10, …
‘StepTolerance’, 1e-8, …
‘FiniteDifferenceStepSize’, 1e-3, …
‘UseParallel’, true);
% Generate a better initial guess
T0 = ones(3, num_projections);
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
for i = 1:num_projections
objFun_i = @(T) computeProjectionError(T, xp(i), yp(i), theta_rad(i), SAD, SID, 1);
estimated_T(:, i) = lsqnonlin(objFun_i, T0(:, i), lb, ub, options);
end
%% Step 6: Calculate RMSE value
% Estimated Tumor Position in X, Y, and Z Direction
Estimate_3D_X = estimated_T(1,:);
Estimate_3D_Y = estimated_T(2,:);
Estimate_3D_Z = estimated_T(3,:);
RMSE_x = Calculate_RMSR (xt,Estimate_3D_X);
RMSE_y = Calculate_RMSR (yt,Estimate_3D_Y);
RMSE_z = Calculate_RMSR (zt,Estimate_3D_Z);
fprintf([‘RMSE for X, Y, Z Direction(mm): ‘ …
‘[%.2f, %.2f, %.2f] mmn’], RMSE_x, RMSE_y, RMSE_z);
%% Step 7: Plot Estimated vs Real Data
Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
toc;
%% Objective function: Computes error between observed and estimated projections
function error = computeProjectionError(T, xp, yp, theta_rad, SAD, SID, num_projections)
for i=1:num_projections
x_est = T(1,:);
y_est = T(2,:);
z_est = T(3,:);
% Compute estimated projection using the perspective transformation
f_theta = (SAD – (x_est(i) .* sin(theta_rad(i)) + z_est(i) .* cos(theta_rad(i))))./ SID;
xp_est(i) = (x_est(i) .* cos(theta_rad(i)) – z_est(i) .* sin(theta_rad(i))) ./ f_theta;
yp_est(i) = y_est(i) ./ f_theta;
end
%% Compute residual error (Method 1)
% Compute the residual error
error_x = xp – xp_est;
error_y = yp – yp_est;
% Return residuals for least squares minimization
error = [error_x; error_y];
end
function RMSE = Calculate_RMSR (True_value,Estimated_value)
% Calculate RMSE
e1 = (True_value(:)-Estimated_value(:)).^2;
RMSE = sqrt(mean(e1,’omitnan’));
end
function Plot_results (Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z)
figure(‘WindowState’, ‘maximized’)
subplot(3,1,1);
plot(Time, Estimate_3D_X,’r’);
hold on;
plot(Time, xt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘X Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,2);
plot(Time, Estimate_3D_Y,’r’);
hold on;
plot(Time, yt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Y Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
subplot(3,1,3);
plot(Time, Estimate_3D_Z,’r’);
hold on;
plot(Time, zt,’b’);
legend (‘Estimate’, ‘Real’);
title (‘Z Dierction’)
xlabel (‘Time (s)’);
ylabel (‘Motion (mm)’)
hold off;
end optimization, fmincon, nonlinear least squares, objective function MATLAB Answers — New Questions
