Seeking Alternative Optimization Methods for a 9-Parameter Optimization Task
Show older comments
Hi there,
I am working on an optimization task where I need to determine 9 optimal parameters that satisfy my objective function (Objective.m). I have experimented with various optimization methods, including fmincon, fminsearch, particleswarm, lsqnonlin, etc, and documented a summary of their performance in computeError_1.m.
So far, fmincon and fminsearch have shown the most promising results, but they are not achieving the level of optimization I need (computeError.m). I am looking for alternative approaches that might improve performance, efficiency, and robustness.
If you have any recommendations or insights on more effective optimization techniques—whether it's hybrid methods, global optimization strategies, or derivative-free approaches—I would greatly appreciate your input.
Thank you in advance for your help!
5 Comments
Did you perform the usual test to see whether your code works as expected: prescribing 9 parameters, producing artificial (slightly disturbed) data with these 9 parameters and then using the respective solver to reproduce the 9 parameters ? Could you include this test case so that we can experiment with it ?
payam samadi
on 31 Mar 2025
William Rose
on 31 Mar 2025
I like the recommendations of @Torsten very much. You suggest, in your reply to him, that you have done what he recommends. However, it seems to me that Test_8.m does not really do what @Torsten recommends. He recommends that the initial guess for the unknown vector be the correct value, or a slightly perturbed value. Test_8.m calls computeError.m, which estimates the 9 unknown parameters with fminsearch (i.e. with simplex algorithm, unconstrained) and with fmincon (constrained). fminsearch() starts from initial guess LvecInit, which is not the actual value of [mx,my,mz,sx,sy,sz,rxy,ryz,rzx] for this data set (Sample_2.xlsx). fmincon() starts from
LvecStart = LvecInit + 0.1 * randn(size(LvecInit));
which is not a slightly perturbed version of the actual vector of parameters for this data set.
Test_8.m does not display the true value of [mx,my,mz,sx,sy,sz,rxy,ryz,rzx] for the data set and does not display the value of [mx,my,mz,sx,sy,sz,rxy,ryz,rzx] found by minimization. I suggest you add that to the code.
Test_8.m does display the covariance matrix of the actual x,y,z data and the covariance matrix of x,y,z data which is computed using the best-fit values for [mx,my,mz,sx,sy,sz,rxy,ryz,rzx]. The comparison of covariance matrices shows that estimated Var(x) is too large by a factor of about 40. The covariance matrices also show the wrong signs for estimated cov(x,y) and cov(z,x). This indicates that the estimated values for rho(x,y) and rho(z,x) have the wrong signs. These are significant errors in estimation of the unknown parameters.
If you modify Test_8.m so that it does what @Torsten suggests, then it will be interesting to see how close it gets to the correct values for the unknown vector [mx,my,mz,sx,sy,sz,rxy,ryz,rzx].
payam samadi
on 1 Apr 2025
Accepted Answer
William Rose
on 31 Mar 2025
Edited: William Rose
on 31 Mar 2025
0 votes
[Edit: Some of the equations in my answer do not look right. I pasted them in. The equations looked fine in the Matlab Answers editing window, but they do not look fine once my answer is posted. Therefore I am editing them. I have not changed the content.]
In your 17 March 2025 posting on Matlab Answers, you included data for a target that moves in a deterministic way in 3 dimensions, as the scanner rotates around the target. The five test data files you attached abouve also show deterministic motion in 3D. The model with 9 adjustable parameters, which you are fitting in this case, is from Poulsen et al., 2009. That model is for a target that moves randomly. The model is that the target is decribed by a 3-D Gaussian distibution. The 3D Gaussian distribution can be thought of as an ellipse in 3D. The 9 parameters are the centroid of the ellipse, the SDs along x,y,z, and the correlations 
, 
, 
. (From these parameters, one can infer the lengths along the principal axes of the ellipse, and the orientation of the 3D ellipse.) The difference between random and deterministric motion is important. If you are using simulated or real data in which the target moves deterministically as the scanner rotates about the target, then the diferent 2D views of the target will not be random samples of the distribution, and one is likely to obtain incorrect estimates for the parameters.
, , . (From these parameters, one can infer the lengths along the principal axes of the ellipse, and the orientation of the 3D ellipse.) The difference between random and deterministric motion is important. If you are using simulated or real data in which the target moves deterministically as the scanner rotates about the target, then the diferent 2D views of the target will not be random samples of the distribution, and one is likely to obtain incorrect estimates for the parameters. Perhaps the script fails to find the correct solution because it gets stuck in a local minimum. You could try different starting points. You could use the mean location, found by a separate method that may use singular value decomposition method, as the initial guess for the first three elements of the 9-dimensional minimization. This is the method used by Poulsen at al., 2008, pp. 1588-1589.
Other ideas for improvement:
- More restrictive upper and lower bounds for the mean y-coordinate, my. ub(my) = f*max(Yp); lb(my) = f*min(Yp); using all projection angles, using all projection angles. f=approximate inverse magnification=SAD/SID.
- More restrictive upper bound and better initial guess for sy. ub(sy)=(ub(my)- lb(my))/2. sy initial guess = f*std(Yp), using all Yp observations. f=approximate inverse magnification=SAD/SID.
- More restrictive upper and lower bounds for mx and mz. ub(mx) = ub(mz) = f*max(abs(Xp)); lb(mx) = lb(mz) = -ub(mx).
- Better initial guess for sx and sz. Let the initial guesses be sx0 and sz0. Three options for the initial guesses are:
- sx0=sz0=s.d.(Xp), using all projections. This may be a considerable overestimate of sx and sz, for a target whose centroid is not at the origin, because, even if the target does not move in 3D at all (i.e. true sx=sz=0), Xp will vary as the source and detector rotate.
- For each projection angle ai, compute the deviation of the projected point from the projection of the centroid, corrected for approximate magnification:
where f=inverse magnification=SAD/SDD and 
=projection of the initial 3D estimate of the centroid. Then
- Use Xp values from projections with angles between ±30° and 180°±30° to estimate sx0, taking magnification into account:
, using only projections for which 
or 
, and where the Xp values are negated when . Use Xp values from projections with angles between 90±30° and 270°±30° to estimate sz0, taking magnification into account: 
, using only projections for which 
or 
, and where the Xp values are negated when . This method is unlikely to yield useful results if there are few projections in each set. - Instead of using initial guesses for sx, sy, sz computed with a specific strategy, as described immediately above, use values chosen at random within the 3D cube defined by the lower and upper bounds for sx, sy,sz. Choose the best fit values that minimize the objective function among all the different initial guesses. This strategy will reduce the chance of getting stuck in a local minimum that is not the global minimum. To sample the different parts of this cube, we probably should use 2^3=8 to 3^3=27 different initial guesses.
- Instead of using initial guess = 0 for rx, ry, rz, use values chosen at random in the 3D cube with bounds ±1. Choose the best fit values that minimize the objective function among all the different initial guesses. Use 8 to 27 different random initial guesses to explore the cube.
- If we combine random initial guesses for sx, sy,sz and for rx, ry, rz, we get a 6D hypercube. To sample the different parts of this hypercube, we probably want 2^6=64 to 3^6=729 different initial guesses. This could take a long time.
22 Comments
William Rose
on 31 Mar 2025
@payam samadi, Here is a little script that plots x,y,z for each of the five test data files you provided. It also computes the values of the nine parameters (3 means, 3 standard deviations, 3 correlations) for each file. The plots show that x,y,z are not varying randomly. The results from running the script are shown below.
>> checkTestData
File N mx my mz sx sy sz rxy ryz rzx
1 1434 +1.51 -7.36 +2.39 2.01 5.66 2.82 -0.94 -0.97 +0.95
2 1448 +0.87 -3.79 +1.67 0.65 5.03 2.33 +0.94 -0.99 -0.92
3 1465 -0.91 -5.92 +0.44 1.67 8.55 1.87 +0.96 -0.94 -0.89
4 1429 -0.18 -4.68 +1.65 0.61 8.20 4.15 -0.79 -0.99 +0.83
5 595 -0.42 -1.42 -1.83 0.71 2.02 1.87 -0.37 -0.61 +0.83
payam samadi
on 1 Apr 2025
Edited: payam samadi
on 1 Apr 2025
Many thanks for you time and comments.
The five test data files you attached abouve also show deterministic motion in 3D. The model with 9 adjustable parameters, which you are fitting in this case, is from Poulsen et al., 2009. That model is for a target that moves randomly.
Ans: Correct. To do this, I project this 3D motion into the 2D space of the imager to obtain the 2D motion observed by the imager.
From these parameters, one can infer the lengths along the principal axes of the ellipse, and the orientation of the 3D ellipse.) The difference between random and deterministric motion is important. If you are using simulated or real data in which the target moves deterministically as the scanner rotates about the target, then the diferent 2D views of the target will not be random samples of the distribution, and one is likely to obtain incorrect estimates for the parameters.
Ans: I am working with real 3D motion data, and I also suspect that projection inconsistencies could be affecting the results. Since I am following Poulsen et al., 2009, I attempted to construct a projection approach suited for the 3D data.
In doing so, I came across an article from the Poulsen et al., 2009 (Appendix) group [1], which appears to use the same dataset as another publication (1, 2). Given this, I believe the projection method (Xp, Yp) described in that article could be applied here.
Another solution would be to consider my own image geometry for p (projection point) and f (projection point), and then modify the equations (A-4 and A-5) on the Objective.m. What do you think?
Additionally, Poulsen et al., 2008 (pp. 1588-1589) appears to simplify the problem by assuming prior knowledge of the mean values in the X, Y, and Z directions. However, this assumption does not fully reflect the true simulation conditions in my case.
Also thanks for the other ideas for improvement. I will go through them.
Ref:
1. Three-dimensional prostate position estimation with a single x-ray imager utilizing the spatial probability density
2. Real-time prostate trajectory estimation with a single imager in arc radiotherapy: a simulation study
William Rose
on 1 Apr 2025
I do not think you need to try different projection equations to convert the 3D position to a 2D projection. I think thsat because, in our previous discussions, I derived one set of equations for the propjection, and uyou used a different set of equations. I compared the results carefully for a set of sample 3D positions, and the different approaches gave exactly the same answers. Therefore I am confident that the equaitons you are using for projections are correct, if you have not changed them since our earlier discussion.
payam samadi
on 2 Apr 2025
To me, it looks like the cone beam magnification factor should be,
SDD/(SAD + z_|| )
not,
SDD/(SAD - z_|| )
That's assuming the z-axis points from the source toward the detector, instead of away from it.
payam samadi
on 2 Apr 2025
Edited: Walter Roberson
on 3 Apr 2025
payam samadi
on 2 Apr 2025
Edited: payam samadi
on 2 Apr 2025
William Rose
on 2 Apr 2025
I have renamed your current file Objecitve.m to Objective20250402.m for my purposes, to avoid confusion with earlier files of the same name. I see that your current Objective.m is very similar to Objective0.m, which I posted on Matlab Answeers on March 11 - except you have removed my name from the comments at the top.
Your latest Objective.m does not use the covariance matrix A, which is computed from the SDs and rho values passed to the function in vector x. My Objective0.m computes B=inv(A). Your script computes B by generating random numbers and using them to make a positive definite matrix.(* See paragraph below.) Your Objective.m does not use the sigmas or the rhos which are included in vector x, which is passed to Objective.m. Therefore the objective function is not utilizing the parameters which the outer minimization routine is trying to optimize. Fix your Objective.m, or use Objective0.m, which I uploaded to this site on March 11.
* You say in the code comments that you generate B the way you do "because the Construct A matrix (Upper Matrix) results in a negative B matrix". Do you mean that B is not necessarily positive deifnite? I would be very suprised if B is not positive definite, as long as the values passed in vector x are reasonable (i.e. as long as the standard deviations are positive and the rho values are between -1 and +1). If A is a valid covariance matrix, then A will be positive deifnite, and if A is positive definite, then B=inv(A) is also positive definite (see here).
Other notes related to your latest Objective.m:
In your current Objective.m, you compute Mahalanobis distance instead of K. In your comments you call this "A.10", but this is only part of eqn A.10 in Poulsen et al 2009. The calculation of K, which I do in Objective0.m, follows Poulsen's A.10. You then incorporate the det(B) term into FMLE in a different way, and your comment says "Compute log-determinant more robustly". My Objective0.m follows the published equations of Poulsen et al. 2009 more closely than your current Objective.m. I have not done the math or other checks to confirm that your "method 2" gives the same answer as my way, but I recommend that you do those checks, if you have doubts.
If your "method 2" (with Mahalanobis distance) does in fact replicate the formulas of Poulsen et al 2009, then I suspect the problem is due to the fact that your data are distributed quite differently from the 3D Gaussian statistics which are the basis of Poulsen's approach. I said this in my March 31 comment above, but I was not as clear as I could have been. You could test my hypothesis by generating simulated 3D data that DOES have a 3D Gaussian distribution. Then you could see if the minimization algorithm which uses your Objective.m (after you fix the problem I identified above) will stop, if the intial guess for vector x matches the parameters of the simulated data.
William Rose
on 3 Apr 2025
The figures in your files Sample_1.xlsx to Sample_5.xlsx show units of mm in the axis labels for x,y,z. You say in the code comments that SAD and SID are 100 and 150 cm. SAD and SID must be expressed the the same units as x,y,z in order to get correct values for projected Xp,Yp. The minimization should still work if the units are mismatched, because the incorrect Xp,Yp values will be mapped back to x,y,z using the same incorrect assumption, theus cancelling out the error. But it would be a good idea to have consistent units, so that the Xp,Yp values are correct.
I think the marker position units in the sample files really are mm, because motions on the order of 10-15 cm would be unrealistically large.
I think the values of 100 and 150 for SAD and SID really are cm, as the code says, because 100 and 150 mm would be too close: the source and detector would be inside the body.
payam samadi
on 3 Apr 2025
Edited: payam samadi
on 3 Apr 2025
Sam Chak
on 3 Apr 2025
It should not take long to locate the file 'Objective0.m' in the recently posted questions. Furthermore, logical deduction suggests that there is a high chance that the original filename was edited and its contents slightly modified. Otherwise, the author’s name would not have been removed from the introductory comments at the top of the function file.
% W. Rose % 20250311
payam samadi
on 3 Apr 2025
Let Y be a Nx3 array of measured 3D positions. The eigenvectors of the covariance matrix of Y are the principal directions of Y. In other words, the eigenvectors (which are guaranteed to be orthogonal) correspond to the axes of the best-fit (in a least-squares sense) ellipsoid that fits the data. The square roots of the eigenvalues are the standard deviations of the data with respect to the principal axes. You can use the array of eigenvectors to rotate the data to an orientation where the principal axes are aligned with the x,y,z axes. (An alternative interpretation is that you can you the array of eigenvectors to express the original data in terms of the rotated coordinate frame defined by the principal axes.) When I applied these ideas to your data, I got a new understanding of your data: I learned that the largest eigenvalue is a lot bigger than the other two, and the smaller eigenvalues are similar. That means the marker moves mainly along one dimension, without much deviation in the orthogonal directions.
The script below applies these ideas to Sample 2.xlsx, which is the data set used by your Test8.m.
Y=readmatrix('Sample 2.xlsx');
Y=Y(:,2:4); % discard column 1
covY=cov(Y); % covariance matrix
[V,D]=eig(covY,"vector"); % eigenvectors and eigenvalues
s=D.^0.5; % square roots of eigenvalues
fprintf('Cov(Y)=\n'); disp(covY)
fprintf('Standard deviations along prinicpal axes=\n'); disp(s')
Display the data:
figure;
subplot(2,2,1); plot(Y(:,1),Y(:,2),'r.')
grid on; xlabel('X'); ylabel('Y'); axis equal
xl=xlim; % get x-axis limits
subplot(2,2,2); plot(Y(:,3),Y(:,2),'r.')
grid on; xlabel('Z'); ylabel('Y'); axis equal
subplot(2,2,3); plot(Y(:,1),Y(:,3),'r.')
grid on; xlabel('X'); ylabel('Z'); axis equal
xlim(xl); % set x-axis limits to match plot above it
The fact that this data varies primarily along one direction may have implications for fitting. For example, you might decide that you only need to find the mean location (3 parameters) and the principal direction of motion (2 parameters) and the standard deviaiton along the principal direction (1 parameter). In this case, you are estimating 6 parameters instead of the 9 parameters in Poulsen's model. When you eliminate poorly constrained parameters from a model, you often find that you can get more reliable estimates for the remaining parameters.
payam samadi
on 4 Apr 2025
Moved: Torsten
on 4 Apr 2025
William Rose
on 4 Apr 2025
Edited: William Rose
on 4 Apr 2025
[edit: fix spelling error]
Thank you for noting that Poulsen 2008 set rLR-CC and rLR-Ap to zero. That is a good example of reducing the number of adjustable parmeters in the model in order to get a better fit for the remaining parameters. I understand that that decision was reasonable for the prostate, but might not be reasonable for a thoracic or abdominal tumor. Thank you also for explaining that Sample 1.xlsx through Sample 5.xlsx are from a liver tumor.
I am impressed that you tried 23 different methods to assure that B is positive definite. But I don't understand how B was not positive definite. How did you determine that B was not positive definite? If A is a covariance matrix, then A and B=inv(A) will be positive definite. Can you give an example of a case when B is not positive definite?
I think I am at risk of forgetting the big picture. Please tell me if the following paragraph is a reasonable sumary of the big picture.
You want to estimate 3D tumor marker location and variation from 2D projections. Fortunately, you have some actual 3D data (Sample [1..5].xlsx), so you can simulate the 2D projections, run your estimation algorithm, and compare the fitted estimates for mean position, stadard deviaitons, and correlations to the actualvalues drived from the 3D data directly. You are not happy with the results you are getting, using the algorithm of Poulsen et al. 2009.
William Rose
on 5 Apr 2025
I implemented the suggestions I made on 17 March. The script successfully estimates the nine marker position and movement parameters, for all five sample files you provided. A screenshot of results is below. The screenshot shows results for each of the five sample data files. For each file, line 1 is the true values of means, SDs, and correlations, derived from the x,y,z data in the sample data files. Line 2 are the estimated values, based on the 2D projections. The agreement is very good, for all 9 parameters, for all five files. The screenshot also shows the distance beween the true mean in 3D and the esitmated mean position. The distance is < 0.2 mm in all cases, which is probably smaller than the marker itself.
The script calls computeError_e(), which is similar to computeError_d. The objective function is Objective0.m (posted here 11 March). Function computeError_e() uses an initial guess and upper and lower limits that are determined from the Xp,Yp data. Specifically, computeError_e() calls computeError_d() and Objective2() to get an initial estimate for the mean position. I emailed Objective2.m to you on 17 March; I see you posted it here, almost identical (but without my name) on 21 March. Function computeError_e() uses the standard deviation of the Xp values, corrected for image magnification, as the initial guess for sd(X) and sd(Z), and it uses the standard deviation of the Yp values, corrected for image magnification, as the initial guess for sd(Y) (as I suggested in my 17 March email and on this site on 31 March). The upper and lower limits for the mean are initial mean estimate +- 2*initial SD estimate. The lower and upper limits for the SDs are 0 and 3*initial SD estimate.
payam samadi
on 5 Apr 2025
William Rose
on 5 Apr 2025
@payam samadi, No problem. Here is the script and functions to get the results above. You have the data files.
I also made 5 simulated data files (attached) that have the same covariance matrices as your 5 data files, but the simulated data has a 3D Gaussian distribution, unlike the real data files which you provided. By uncommenting one line, and commenting out the other line, in the main script, the script will use the simulated files instead of the real data files. If you do this, you will see that Lestim (the 1x9 vector estimated from 2D projections) is even closer to Ltrue (the vector computed from the X,Y,Z data) than when we analyze real data. I suspect the estimation is better with the simulated data because Poulsen's method is a maximum likelihood estimator for Gaussian data. The simulated data IS Gaussian, and the real data is not. But, even with real data, Lestim is very close to Ltrue.
Thank you for your kind comments. I too have learned a lot on this site from the others who know more about math and Matlab than me.
payam samadi
on 6 Apr 2025
Edited: payam samadi
on 6 Apr 2025
@Payam,
The 1x9 vector L, which we estimate from 2D projection data using the scripts I provided above, contains information about the location and motion of the target. Here is a script that uses the L vector (Lestim1, estimated from 2D projections of Sample 1.xlsx, see previously posted code) to estimate the principal axis of movement and the percent of marker location variation that is accounted for by movement along the principal axis. I assume that +X=left, +Y=cranial, +Z=anterior.
visualize3dGaussian0
payam samadi
on 6 Apr 2025
Here is a revised version of the script to assist in visualizing the 3D Gaussian distribution. The input is the 1x9 L vector = [mean(x),mean(y),mean(z),sd(x),sd(y),sd(z),corr(x,y),corr(y,z),corr(z,x)].
L may be estimated from 2D projections or may be computed directly, if the x,y,z coordinates are known.
The script displays the orientation of the principal axis of the 3D Gaussian distribution, and displays the fraction of total variance accounted for by movement along the principal axis.
The script generates 2 figures showing an ellipsoid that bounds 90% of the probability for the marker location. Change the value 0.9 in the code, to generate an ellipse that bounds a different probability.
The first figure is designed to be rotated as desired by the user. The second figure shows three orthogonal two-D views of the ellipsoid, for a supine patient. The ellipsoid color has no specific meaning, but it helps the viewer know how the different views are related.
visualize3dGaussian1
More Answers (0)
Categories
Find more on Get Started with Curve Fitting Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
