set of linear equations with uncertainty
in [Matlab]
|
From: David on 10 Aug 2010 10:45 Greg Heath <heath(a)alumni.brown.edu> wrote in message <871417f8-d089-40f5-ae39-0922bd514fab(a)f42g2000yqn.googlegroups.com>... > On Aug 9, 4:19 pm, "David " <dmkl...(a)mdanderson.org> wrote: > > Hello, > > > > I have two measurement matrices - A and b - with which I can set up the matrix equation Ax = b. I have been solving for the x matrix in MATLAB by using the following: > > > > x = pinv(A)*b, > > > > and so far this has given me pretty good results. However, thus far I have not taken the measurement uncertainties in A and b into account. I admit that linear algebra is not my strong point, and would like to ask the Newsgroup for suggestions on solving the following using MATLAB: > > > > (A + δA)(x + δx) = (b + δb), > > > > for x and δx, where δA. δx, and δb are the uncertainty matrices of A, x, and b, respectively. > > > > Any help would be greatly appreciated, and thanks ahead of time! > > Dave > > Please change your font so that sequences of weird parameters > are not posted on the newgroup. > > About 50 years ago I had a homework problem that required > the following proof for invertible A: > > ||dx||/||x|| <= cond(A)*(||dA||/||A|| + ||db||/||b||) > > This can be derived by linearizing and combining the equations > > 1. A*x = b > 2.||b|| <= ||A||*||x|| > 3. dx = inv(A)*( db -dA*x) > > Hope this helps. > > Greg Hi Greg, Sorry for the font problems, the deltas look fine on my machine. I should've just used dA, etc. Yeah, I've calculated an upper bound for the relative error in the x matrix using the condition number and matrix norms. I'd just like to see if I could get uncertainty values corresponding to each term in x. Thanks, Dave
From: David on 10 Aug 2010 14:32 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <i3r0ml$92s$1(a)fred.mathworks.com>... > > > > > Hi Bruno, Thanks for your fast response. > > > > To give a little more detail, both A and b are m × 2 matrices, and the x matrix is 2 × 1. > > The dimensions do not match, > - if A is m x 2, and b is m x 2, then x must be 2 x 2. > - if A is m x 2, and x is 2 x 1, then b must be m x 1. > > >I usually use m = 2, but would like to investigate the effects of using A and b having m > 2. The error in the b matrix is indeed Gaussian. The distribution of the uncertainty in A is technically unknown, but practically I think it can be assumed to be Gaussian without causing any problems. I'd ultimately like to get an uncertainty term for each value of the x matrix. > > You still don't fully answer my question, which is: how do you input the uncertainty of A and b? The Gaussian is characterized by the covariance *matrix*, which is a matrix of the dimension of the vectors. Do you know them? > > The simplest case is the covariance matrix is the fixed variance times identity, in this case use total least-squares like this (assuming b is a column vector): > > % Data > A = rand(3) > b = rand(3,1) > > % standard deviation of A, the scalar here is simple case, > % this can get more complicated: a 9 x 9 sdp matrix (6DOFs) needed to be entered > sigmaA= 0.1 > % standard deviation of b, similar remark to above > sigmab = 1 > > % engine > M = [1/sigmaA*A 1/sigmab*b] > [U S V] = svd(M,0); > n = size(A,2) > x = -V(1:n,end)/V(end) > > % Bruno Hello again Bruno, I mistyped in my first reply to you. The b matrix is m x 1. With respect to the form of the uncertainty, I have uncertainty terms corresponding to each of the terms in A and b. I'm pretty rusty when it comes to the covariance matrix, but I'm looking into it. Each of the terms in A and b are averages, and each of the terms in dA and db are the standard deviations of the data that make up those averages. For example, let's say that A = [mean([s1,s2,s3]), mean([t1,t2,t3]); mean([s4,s5,s6]), mean([t4,t5,t6]); ...] then dA = [std([s1,s2,s3]), std([t1,t2,t3]); std([s4,s5,s6]), std([t4,t5,t6]); ...] (Actually, the terms in A result from median calculations, and the terms in dA result from an estimation of the uncertainty in those median calculations, but for this discussion I guess we can just use mean for A and std for dA.) Likewise, matrix b = [mean([u1,u2,u3]); mean([u4,u5,u6]);...] and db = [std([u1,u2,u3]); std([u4,u5,u6]);...] where s, t and u are just dummy scalars. Using these I'd like to obtain a solution matrix x and a identically-sized dx matrix in which each term is the uncertainty of the corresponding term in x when m >=2. In the end I will be using the results to solve the equation s*x1 + t*x2 = u, where s and t are measured values, u is unknown, and x1 and x2 are the coefficients determined using measurements corresponding to reference u values, as in a calibration procedure. My goal is to account for as many sources of error/uncertainty as I can. Intuitively, I would think that the terms in dx would decrease with increasing m (assuming the relation is truly linear), due to the presence of more data to apply the least-squares fit to. I hope that's a bit more clear. Thanks again for your help, Dave
From: Bruno Luong on 11 Aug 2010 03:05 Well David, you ask indeed two separate questions here, and I'm not sure which one you need in priority. 1. Given the uncertainty of matrix/data, what is the right way to take into account them and solve A*x = b? The answer is WTLS. WTLS in general does not have close-form solution. However it can be solved by numerical optimization. In some simple cases, such as the one I wrote earlier, when the uncertainties are simple scalar times identity, the WTLS can then be solved using SVD, similar to TLS. In your case I'm not sure you can use SVD. I believe the assumption that makes WTLS transformable to TLS is that the covariance matrix of A can be written as tensorial-form corresponding to row/column, which so far you haven't stated (?). You just need to work out from literature of WTLS (see also Makkovsky / Huffel beside the references I gave earlier), but I strongly suggest to acquire a more solid knowledge on covariance matrix prior to tackle the reading task. 2. How to estimate the error of x (Greg gives a bound using condition number, but I believe it's not accurate by any mean)? Again, you seem to think that the error is the same characterize by a vector of the same size (by saying dx is the same size as x), which is not true in general. The error is characterized by a covariance *matrix*. The "error bar" of the estimated parameters is square-root of the diagonal of this matrix, and it only tells half of the story at best. Furthermore the TLS and WTLS is non-linear estimator (contrary to LSQ), hence the error dx does not even follow Gaussian distribution. The covariance matrix is hard to estimate. People usually compute an approximation of it. For LSQ inversion, the computation is already complicated when take into account dA. The Jacobian of x wrt A is given by Golub/Pereyra formula. I believe there is an approximate formula for TLS given by Schaffrin 2006. I'm not aware of paper that give a practical formula for WTLS. As you see, the problem you state is kind of hard to solve generally and the field is still being developed as we speak. I suggest you to state more clearly what you want (1 or 2?), and may be relax a little bit what you want to solve, because it can get nasty if you don't fix a limit. So far it seems that you are already satisfied with LSQ, which corresponds to WTLS with dA = 0 and db = cte2*eye(length(b)). Without going to a generic framework, you should get better result by simply relaxing dA = 0 to dA = cte1*eye(numel(A)). Hope it's not confusing. Bruno
From: David on 11 Aug 2010 10:38
Hi Bruno, Thank you very much for your thorough reply. I guess I am dealing with two different problems. For now I think I'll make do with the upper bound for relative error that I get using the condition number - it's around 1%, which is not too bad for my application. However, I will continue looking through the literature. Thanks for all of the references - I'm looking through one of Chris Paige's papers now. Cheers! Dave |