Optimization method based on the Nonlinear least squares as well as fmincon for the defined objective function (Non-linear function).

10 visualizaciones (últimos 30 días)
payam samadi
payam samadi el 17 de Mzo. de 2025
Editada: Walter Roberson el 29 de Mzo. de 2025
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] mm\n'], 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
  3 comentarios
Mostrar 1 comentario más antiguo
Matt J
Matt J el 17 de Mzo. de 2025
Since you have no constraints, you may as well use fminunc instead of fmincon, or fsolve instead of lsqnonlin.
payam samadi
payam samadi el 18 de Mzo. de 2025
I made a change in step 5 and defined a new function called "computeProjectionError_b." However, I got some unexpected results that are quite confusing. I'm not sure where I went wrong, but here's what I'm trying to do: In the first loop, I used only the initial guess. For the following iterations, I used the estimated results to update the variables. I should mention that this might not be the best approach, but it was the idea that came to mind. Any help or alternative methods would be greatly appreciated.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuesta aceptada

Matt J
Matt J el 18 de Mzo. de 2025
Editada: Matt J el 18 de Mzo. de 2025
Ran in:
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.
But for every time Time(i), you appear to have only one projection view of the tumor position. This is not enough, I'm afraid. As noted by Torsten as well, it will be an underdetermined problem. However, since this appears to be a radiotherapy scenario, and since radiotherapy treatment systems typically have two x-ray imagers, you should be able to get two projections per tumor position. If we adapt your code to this scenario (see below), we can see the inverse problem is easily solved:
%% 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=height(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)';
theta_rad = deg2rad(alpha)+[0,pi/2];
num_detectors=width(theta_rad);
%% Step 3: Projection Model: Compute 2D projections
[Xp,Yp]=deal(nan(num_projections,num_detectors)); % allocate arrays
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (true_T(i,1) .* sin(theta_rad(i,j)) + true_T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
Xp(i,j) = (true_T(i,1) .* cos(theta_rad(i,j)) - true_T(i,3) .* sin(theta_rad(i,j))) ./ f_theta;
Yp(i,j) = true_T(i,2) ./ f_theta;
end
end
xp=Xp; yp=Yp;
%% Step 4: Define the objective function
tic
T = optimvar('T',[numT,3]);
[eqn1,eqn2]=deal(cell(num_projections,1));
for i=1:numT
for j=1:num_detectors
f_theta = (SAD - (T(i,1) .* sin(theta_rad(i,j)) + T(i,3) .* cos(theta_rad(i,j)))) ./ SID;
eqn1{i,j}=xp(i,j).*f_theta == (T(i,1) .* cos(theta_rad(i,j)) - T(i,3) .* sin(theta_rad(i,j))) ;
eqn2{i,j}=yp(i,j).*f_theta == T(i,2) ;
end
end
prob=eqnproblem();
prob.Equations.eqn1=vertcat(eqn1{:});
prob.Equations.eqn2=vertcat(eqn2{:});
s=prob2struct(prob);
Linear equation problem is overdetermined. Equations will be solved in a least squares sense.
estimated_T=reshape(s.C\s.d,[],3);
toc
Elapsed time is 5.777380 seconds.
%% 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] mm\n'], RMSE_x, RMSE_y, RMSE_z);
RMSE for X, Y, Z Direction(mm): [0.00, 0.00, 0.00] mm
%% Step 7: Plot Estimated vs Real Data
Plot_results_Matt(Time, xt, yt, zt, Estimate_3D_X, Estimate_3D_Y, Estimate_3D_Z);
●
  34 comentarios
Mostrar 32 comentarios más antiguos
Matt J
Matt J el 24 de Mzo. de 2025
Editada: Matt J el 24 de Mzo. de 2025
Do you think other power optimization methods such as GA and gravity search could handle and optimize these nine parameters?
It might help with the local minima, but I feel like this problem should be easy to come up with a good Lvecinit. I don't know how you are deriving your Lvecinit currently.
Still, the results are so far from our expectations.
In Objective_2 you have imposed bounds on x(4:9) which I said you should not have. You can impose bounds on x(1:3), but I don't know where you are getting them. In your test data, it looks like the target point drifts more than a centimeter from isocenter, so +/-10 is the least I would choose.
Aside from all that, you still don't appear to have run the experiment initializing at ground truth. I will take a break from this thread until you do.
payam samadi
payam samadi el 25 de Mzo. de 2025
Editada: Walter Roberson el 29 de Mzo. de 2025
Dear Matt,
Many thanks for your reply.
To be honest, I got completely confused.
I don't know how you are deriving your Lvecinit currently.
I just considered an initial guess for Lvec. Also, I found this from this article. 
For each CBCT scan, three optimizations with three sets of starting points for (σ LR, σ CC, σ AP) and (ρ LR-CC, ρ LR-AP, and ρ CC-AP) were performed and the one resulting in the highest likelihood was selected.
Therefore, I used the Multi-start strategy (attached computeError.m here): Run multiple times with different initial guesses for Avoiding Local Minima (the following code is added to the computeError.m):
LvecInit=[5,1,0,2,8,1,0,5,-4]; % initial guess for Lvec
ub = 2* ones (1,9); % lower bounds
lb = -2*ones (1,9); % upper bounds
% For Method 4 (Method 4 + fmincon)
options = optimoptions('fmincon', ...
'Algorithm', 'sqp', ... % Algorithm choice 'sqp', 'sqp-legacy', 'interior-point', 'active-set'
'MaxIterations', 1000, ... % Maximum number of iterations
'MaxFunctionEvaluations', 5000, ... % Maximum function evaluations
'OptimalityTolerance', 1e-12, ... % Tolerance for stopping
'StepTolerance', 1e-12, ... % Minimum step size before stopping
'Display', 'iter', ... % Display output at each iteration
'UseParallel', true);
for i=1:N
% Find 1x9 vector Lvec that minimizes F_MLE for marker i.
xpv = Xp(:,i); % column vector, length Np
ypv = Yp(:,i); % column vector, length Np
fun=@(x) Objective_3(x, xpv, ypv, SAD, SID, alpha); % function to
% minimize (Included matrix B is positive definite) + Recommended by Matt + Modified by Payam.
%% Methods to Optimize Objective Function
% Multi-start strategy: Run multiple times with different initial guesses
num_restarts = 5; % Number of restarts
best_Lvec = [];
best_fval = Inf;
for j = 1:num_restarts
% Run Nelder-Mead optimization
% [Lvec, fval] = fminsearch(fun, LvecInit, options);
% [Lvec, fval] = fminunc(fun, LvecInit, options);
if j == 1
[Lvec, fval] = fmincon(fun, LvecInit, [], [], [], [], lb, ub, [], options); % Method 4 + fmincon
else
[Lvec, fval] = fmincon(fun, best_Lvec, [], [], [], [], lb, ub, [], options); % Method 4 + fmincon
end
% Store best solution found
if fval < best_fval
best_fval = fval;
best_Lvec = Lvec;
end
end
LvecArray(i,:) = best_Lvec; % save Lvec in row i of LvecArray
end
For the first iteration, the model runs, but for the next iterations, it gives:
Local minimum possible. Constraints satisfied.
fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are
satisfied to within the value of the constraint tolerance.
My propsoed strategy is correct for this statement.
For each CBCT scan, three optimizations with three sets of starting points for (σ LR, σ CC, σ AP) and (ρ LR-CC, ρ LR-AP, and ρ CC-AP) were performed and the one resulting in the highest likelihood was selected.
Also, I told you that you need to set your lb, ub to bound to ensure that -1<=[rxy,ryz,rzx]<=1.
I applied this to Objective_3.m (attached here) with different methods to make sure that the Matrix B is positive definite (please see Objective_3.m).
mx=x(1); my=x(2); mz=x(3); % extract values from x
sdx=x(4); sdy=x(5); sdz=x(6);
rxy=x(7); ryz=x(8); rzx=x(9);
%% Ensure Proper Bounds on Correlations (ρxy, ρyz, ρzx)
% If the correlation coefficients rxy, ryz, and rzx are outside the
% valid range (-1 ≤ ρ ≤ 1), A may become a non-positive definite:
rxy = max(min(rxy, 1), -1);
ryz = max(min(ryz, 1), -1);
rzx = max(min(rzx, 1), -1);
%% Calculate FMLE
r0=[mx;my;mz]; % theoretical mean position
Np=length(xp); % number of projections
% Construct the covariance matrix A
A=[sdx^2, rxy*sdx*sdy, rzx*sdz*sdx;...
rxy*sdx*sdy, sdy^2, ryz*sdy*sdz;...
rzx*sdz*sdx, ryz*sdy*sdz, sdz^2];
% Method 6: Combining Multiple Techniques
A = (A + A') / 2; % Method 2
[~, p] = chol(A); % Method 5
if p > 0 % If A is not positive definite
A = A + (abs(min(eig(A))) + 1e-1) * eye(3); % Increase diagonal elements
end
A = A + 1e-10 * eye(3); % Step 1: Final regularization
I think, for now, I see the problem with the local minima.

Iniciar sesión para comentar.

Más respuestas (2)

Torsten
Torsten el 17 de Mzo. de 2025
Editada: Torsten el 17 de Mzo. de 2025
You can add a multiple of the vector you obtain by the below command "double(null(A))" to "estimated_T", and you will still get a solution to your linear system of two equations. The system is underdetermined and there exist infinitly many solutions for it.
% @ 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';
% Run optimization for each projection
estimated_T = zeros(3, num_projections);
T = sym('T',[3 1]);
for i = 1:num_projections
f_theta = (SAD - (T(1) .* sin(theta_rad(i)) + T(3) .* cos(theta_rad(i))))./ SID;
eqn1 = xp(i)*f_theta - (T(1) .* cos(theta_rad(i)) - T(3) .* sin(theta_rad(i)))==0;
eqn2 = yp(i)*f_theta - T(2)==0;
[A,b] = equationsToMatrix([eqn1,eqn2]);
double(null(A))
estimated_T(:, i) = double(A)\double(b);
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] mm\n'], 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;
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
  6 comentarios
Mostrar 4 comentarios más antiguos
Torsten
Torsten el 18 de Mzo. de 2025
Editada: Torsten el 18 de Mzo. de 2025
I noticed that I made a mistake and the number "1" should be "num_projections".
After replacing 1 by num_projections, your code throws an error. In function "computeProjectionError", you try to access theta_rad(i) for 1 <= i <= num_projections in your loop, but you only pass theta_rad(i) to the function which is a single scalar value.
payam samadi
payam samadi el 18 de Mzo. de 2025
I made a change in step 5 and defined a new function called "computeProjectionError_b." However, I got some unexpected results that are quite confusing. I'm not sure where I went wrong, but here's what I'm trying to do: In the first loop, I used only the initial guess. For the following iterations, I used the estimated results to update the variables. I should mention that this might not be the best approach, but it was the idea that came to mind.
Any help or alternative methods would be greatly appreciate

Iniciar sesión para comentar.

payam samadi
payam samadi el 28 de Mzo. de 2025
Editada: payam samadi el 28 de Mzo. de 2025
Dear Matt,
I have carefully reviewed the code and made several refinements. While I observed noticeable improvements compared to the previous version, the results are still not fully meeting my expectations. Below is a summary of my findings and ongoing challenges:
1. Ensuring Matrix B is Positive Definite:
  • I tested 24 different methods to enforce positive definiteness on Matrix B. However, only methods 9 and 10 produced meaningful results (see Objective_3.m, lines 57-213).
  • To validate this, I had a colleague test the code using a Gravitational Search Algorithm (GSA) for optimization. Interestingly, modifying the method for ensuring B's positive definiteness led to wildly different results (Test 6-GSA.zip file).
  • Although multiple methods exist (to make sure Matrix B is positive definite), methods 9 and 10 performed best (though not optimally). However, since they rely on L=zeros(3) or L = tril(randn(3)), enforcing proper bounds on correlation coefficients (ρxy, ρyz, ρzx) becomes ineffective.
2. Improved Mahalanobis-Based Approach (Method 2):
  • To address numerical instability due to high condition numbers (nearly singular B), I implemented an improved Mahalanobis-based approach (see Objective_3.m, lines 236-239).
  • This modification aims to stabilize computations and enhance reliability.
  1. Initial Guess for LvecInit:
  • Currently, I am using [0,0,0,1,1,1,-0.9,-0.9,0.9] as the initial guess across different samples [1-5].
  • While this yields reasonable results, further refinement may improve convergence and accuracy.
3. Bounding Mean Tumor Motion (mx, my, mz):
  • I incorporated constraints into Objective_3.m to ensure that mean tumor motion (mx, my, mz) remains within the range -10 ≤ m ≤ 10.
  • This helps maintain realistic motion estimations while improving numerical stability.
4. To avoid getting into the local minima:
  • Two functions computeError and computeError_1 are proposed to avoid getting into the local minima. The computeError_1 is better, and included two different optimzation; strating with fminsearch and then for the next iteration fmincon is used.
Overall, these updates have resulted in improved stability, but there is still room for this code. Let me know your thoughts on the next steps or if you have any suggestions.
Here is results of the code

Categorías

Más información sobre Solver Outputs and Iterative Display en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by