I am currently trying to implement the method described in this paper. In short we have a system of the form $a=B\times c$. Where

$$a_i = \int d^3r \space w(r)f_i(r)t(r)\quad\text{and}\quad B_{ij} = \int d^3r \space w(r)f_i(r)f_j(r).$$

The solution is given by $c = B^{-1}a$. In this case $w,t,f_i$ are volumes, specifically a sequence of images obtained from an x-ray simulation. When I solve these integrals using MATLAB's trapz function I get a poorly-conditioned $B$ matrix.

  • Is there any other way to integrate this type of data (images)?
  • Or a different way to calculate my $a$ and $B$ terms?

You can try to improve one both sides. One issue is that smooth separable kernels $w$ or integration scheme like trapz tend to reduce rank, and pose conditioning problems. Did you use pinv or X = (A)\(b) already?

  1. You can invest on non-separable or higher order integration rules, see for instance n-dimensional simplex quadrature;
  2. You can invest in regularization, for instance with the simple diagonal update: $B_c = B + c*\operatorname{eye}(\operatorname{size}(B))$. If this is not sufficient, there are quantities of constrained or penalized methods (LASSO for instance), implemented for instance in the CVX toolbox (e.g. Premade solvers for specific problems (vector variables)).
| improve this answer | |
  • 1
    $\begingroup$ Yes, there is nothing to be worried about using pseudo inversion in these cases. Maybe the reason for which you get an ill conditioned matrix is not trapz, but maybe the data are not sufficiently informative to yield a single solution. $\endgroup$ – LJSilver Oct 20 '16 at 8:18
  • $\begingroup$ I have tried regularization with very little success. Will now attempt the other suggestion. However, how could my data not be sufficiently informative? $\endgroup$ – RCountZero Oct 20 '16 at 10:58
  • $\begingroup$ Well in the case in which your matrix is not invertible you just have infinite solutions in the least squares sense. Just pick one $\endgroup$ – LJSilver Oct 20 '16 at 22:32
  • 1
    $\begingroup$ Regularization worked. I just had to choose the correct number. Thank you all for your help. $\endgroup$ – RCountZero Oct 21 '16 at 6:51
  • $\begingroup$ Good news, the typical issue with regularization is indeed the "correct number" $\endgroup$ – Laurent Duval Oct 21 '16 at 7:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.