Incorrect inversion of (complex) matrices

[ajlawren]

Hello,
This is an amazing package, I am thoroughly impressed. Thank you for all the time and effort!

I am trying to inv a square matrix A composed of Particle objects:

Using MonteCarloMeasurements
A_bar = randn(5,5)
A_std = abs.(randn(5,5))
A = A_bar .± A_std 

Calculating the inverse of this with MonteCarloMeasurements (with the @unsafe mode) yields

@unsafe pmean.(inv(A))
5×5 Matrix{Float64}:
  0.991553    1.16502     -1.94536   -1.50509    0.439208
  0.232711   -0.00783246   0.414581  -0.0993158  0.409937
  1.1345      0.741083    -1.54862   -0.970953   0.749576
 -0.794642   -0.256104     0.906283   0.706037   0.0634964
  0.0174808   0.666404    -2.28201   -1.68       0.617051

while calculating the inverse of the mean yields

inv(A_bar)
5×5 Matrix{Float64}:
 -0.0441282    0.207341    0.382708   0.243892   -0.160185
 -0.137243    -0.0641063   0.438657  -0.0528539   0.617828
 -0.00204959  -0.114085    0.484008   0.742596    0.079385
  0.178726     0.340945   -0.502247  -0.466079    0.238289
 -0.934811    -0.156226   -0.208956  -0.0631625   0.273484

Is this a result of the @unsafe mode? Or is this impossible with Monte Carlo methods?

This operation is handled fine with Measurements, interestingly enough. I am ultimately interested in using inv on complex uncertain matrices, but this is the more fundamental issue.

Thank you!

[baggepinnen]

Hello and thank you for your kind words!

Computing the mean of the result mean(f(x)) is in general not the same as computing the result of the mean f(mean(x)) which is what Measurement.jl does. These two are only the same for linear functions, but this is not true for nonlinear functions. See this past discussion with a similar example
https://discourse.julialang.org/t/uncertainty-propagation-using-montecarlomeasurements-gives-weird-results/124745

[ajlawren]

Hi @baggepinnen,
Thank you very much for your informative response! I've been thinking about this since you posted it, and yes I see what you mean and that's definitely what I'm experiencing. It's interesting that matrix inversion has this "property" of nonlinear functions, when it's literally a linear computation.

I'm not sure about how to proceed in this instance with MonteCarloMeasurements, when I have a computation that depends on the use of mean(). Does that not fit in this theoretical framework? Do I need to reframe the problem for it to work in MonteCarloMeasurements?

For the record, I'm trying to get the uncertainty of a power spectral density matrix and compute the FRF transfer function from Bendat and Piersol.

[baggepinnen]

Matrix inversion is not a linear function, you can easily convince yourself of this by considering inversion of a scalar, the graph of $1/x$ is not a straight line through the origin.

It's not clear to me what you are trying to do. When performing monte carlo sampling, you draw a number of samples from the distribution of the parameters, perform the same computation on all of them, and then compute some statistic that you are interested in of the results. Which of these three steps do you have problems with? You never compute any summarizing statistic of intermediate results, for example, if the inverse matrix above is not what you are ultimately interested in, don't bother computing the mean of it.