I have a linear system of equations in the form Ax = b. The elements of A and b were experimentally determined and as such have some uncertainty. Within each row, A and b are correlated. Between rows, A and b are uncorrelated. How would I propagate this to the elements of x?

\begin{equation} \begin{bmatrix} a_{11} \pm \sigma_{a_{11}} & a_{12} \pm \sigma_{a_{12}} & a_{13} \pm \sigma_{a_{13}} \\ a_{21} \pm \sigma_{a_{21}} & a_{22} \pm \sigma_{a_{22}} & a_{23} \pm \sigma_{a_{23}} \\ \end{bmatrix} \begin{bmatrix} x_{11} \pm \sigma_{x_{11}} \\ x_{21} \pm \sigma_{x_{21}} \\ x_{31} \pm \sigma_{x_{31}} \\ \end{bmatrix} =\begin{bmatrix} b_{11} \pm \sigma_{b_{11}} \\ b_{21} \pm \sigma_{b_{21}} \end{bmatrix} \end{equation}

, where $x_{11}+x_{21}+x_{31}=1$, and $x_{11}$, $x_{21}$, $x_{21}$ $>0$. I am trying to find $\sigma_{x_{11}}$, $\sigma_{x_{21}}$, and $\sigma_{x_{31}}$.

I have been looking at the threads on propagation of uncertainty through a linear system of equations and variance of dependent variables but I am unable to combine the two ideas.

There has been a similar unanswered question. Will be hugely grateful for any pointers!

Mark Stone has been extremely helpful with a statistical sampling method. If there is a deterministic way, or some understanding of the challenges of a deterministic way, that would also be of interest, possibly not to just myself.

Update: I tested the solution to this problem but it raised a couple of questions on its own.

  • $\begingroup$ Do you have 2 linear equations in 3 unknowns, and therefore, leaving aside the uncertainty, an under-determined system, not over-determined? $\endgroup$ – Mark L. Stone Mar 8 '16 at 1:05
  • $\begingroup$ Thanks Mark for your comment. Apologies, you may be right! It's very late at night here. It is indeed 2 linear equations. I have 3 unknowns, but the third unknown is a linear combination of the first two. Does it qualify as an independent third unknown? $\endgroup$ – Pete Mar 8 '16 at 1:21
  • $\begingroup$ I don't understand what constraint or relation you have on the 3rd unknown, or as some linear combination of the first two. Really need to understand what the problem is first, then can worry about the solution. $\endgroup$ – Mark L. Stone Mar 8 '16 at 1:23
  • $\begingroup$ Thank you for asking! B is an expression of the extent of x-ray CT absorption within a volume, as a function of pure components A. Here, x represents volume fractions of the components of A, and must sum to 1 (thus the constraint on the third unknown). In this problem, x represents fractions of dolomite, calcite and air; 'A' represent the experimentally measured absorptions of the three pure components; and 'B' represent the volume-averaged absorption at two different conditions where the pure components A respond differently. Please let me know if I may further clarify. $\endgroup$ – Pete Mar 8 '16 at 1:33
  • $\begingroup$ Constraint that the x's sum to 1 is your 3rd linear equation, which has no uncertainty in LHS or RHS elements. Leaving aside the uncertainty, your system is exactly determined, presuming your A instantiation (with 3rd row) is not singular.The most straightforward approach, which allows for any form of dependency, is to perform n stochastic simulation replications, where for each replication, you draw random values for A and B, and solve A_augmented x = B_augmented, where augmented means adding the row specifying x's sum to 1.Then you have n sample values of the x vector, w/ all dependencies. $\endgroup$ – Mark L. Stone Mar 8 '16 at 1:49

Here is a stochastic simulation approach. I am not ruling out that someone may come up with an approach not requiring simulation. Nevertheless, the simulation run time should be negligible, so if the goal is a practical solution, this should qualify.

You have a system of 2 linear equations in 3 unknowns, $A x = b$, augmented by the additional linear equation that the sum of the x's equals 1. This adds a 3rd row to A consisting of all 1's, and a 3rd element of b equal to 1. There is no uncertainty associated with this 3rd row. Let $A_{augmented}$ be the 2 rows of A augmented by a (3rd row) row of 1's, and $b_{augmented}$ be the 2 elements of b with a 3rd element equal to 1. So your complete system of linear equations is $A_{augmented} x = b_{augmented}$

The first step is to either

1) fit distributions to the values of A and b, accounting for any correlations among the elements, including across A and b, then for each simulation replication, draw a sample of A and b from the fitted distributions (for instance, Multivariate Normal which can allow correlation between A and b) or

2) use raw measurements of A and b "as is", with each measurement (combination between A and b) being used for a simulation replication.

For each of n simulation replications, solve the equation $A_{augmented} x = b_{augmented}$, which provides the simulated value of x for that simulation replication. Once you have completed the simulation replications, you can assess the error distribution of the n sample values of x (including correlation across elements). The accuracy you achieve in assessing the uncertainty of x will ultimately be limited by the number m of A and b measurements (decreasing as its square root), and if fitting distributions to A and b, the validity of the distributions you fit. Nevertheless, if using fitted distributions for A and b, you might as well do a large number of replications n, so that you get maximum mileage out of the data you do have - no need for number of replications n to be limited to number of observations m.

  • $\begingroup$ Mark, thanks for this, it is working as far as I can tell. However, I'm getting negative values of x. Since x is a physical volume fraction, how can I add a constraint for only positive values of x? $\endgroup$ – Pete Mar 8 '16 at 17:34
  • $\begingroup$ Instead of solving the system of linear equations, you could do bound-constrained linear least squares, imposing the bound constraint x >= 0. You can use linearly-constrained linear least squares, which will include bound constraints as a special case. You'll have to decide if you want to impose the x's summing to 1 as a hard constraint (which would be a linear equality constraint to the linear least squares problem), or if you want to allow some wiggle room on it. If the former, remove it from the least squares part, and impose as a linear constraint, in addition to the x >= 0 constraint. $\endgroup$ – Mark L. Stone Mar 8 '16 at 17:52
  • $\begingroup$ Thanks Mark. I simulated x and using your method, attempted to back-derive x and its error or standard deviation. I get the right mean but the standard deviation differs significantly. I simply don't understand what is wrong. If you have any comments, or pointers, I'd be immensely grateful. The script was written in Matlab. This is similar to another question posted on this site. $\endgroup$ – Pete Mar 9 '16 at 1:20
  • $\begingroup$ I don't know what you're doing on the APPROXIMATE analytical error propagation, but surely you're not accounting for the nonnegativity constraint, likely among other things. $\endgroup$ – Mark L. Stone Mar 9 '16 at 1:42
  • $\begingroup$ Thanks Mark for this. Do you mean that the constraints affect the propagation of uncertainties of A and x to b in $a_{11} x_{11} + a_{12} x_{21} + a_{13} x_{31} = b_{11}$, and $\sigma_{a11x11} = a_{11} x_{11} \sqrt{\left(\frac{\sigma_{a11}}{a_{11}}\right)^2 + \left(\frac{\sigma_{x11}}{x_{11}}\right)^2} etc$, and $\sigma_{b11}=\sqrt{\sigma_{a11x11}^2 + \sigma_{a12x21}^2 + \sigma_{a13x31}^2}$ is not valid? $\endgroup$ – Pete Mar 9 '16 at 2:09

Basically, the question is how to calculate uncertainties of the derived coefficients in errors-in-variables models.

One way is to use Monte Carlo simulations, as showed in Mark's answer. For more analytical solutions, you may need to read https://en.wikipedia.org/wiki/Errors-in-variables_models and other materials. I haven't used this before. Hope other people can comment on it. You may refer to this post for the difference between uncertainty from MC simulations and uncertainty from analytical derivations: Is Monte Carlo uncertainty estimation equivalent to analytical error propagation?

By the way, you said that 'Within each row, A and b are correlated'. Of course, A1 and b1 is are correlated, otherwise we won't derive X, right? By stating 'Between rows, A and b are uncorrelated', do you mean your observations of (A1, b1) and (A2, b2) are uncorrelated?

  • $\begingroup$ Thanks yuqian. I had a look at the wikipedia page, and I admit it was a struggle. I am primarily an experimentalist and so relegated to simpler maths. Yes (A1, b1) and (A2, b2) are uncorrelated, they come from two independent measurements. $\endgroup$ – Pete Mar 9 '16 at 0:22

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.