A universal and fast method to solve linear systems with correlated coefficients using weighted total least squares

Especially in metrology and geodesy, but also in many other disciplines, the solution of overdetermined linear systems of the form Ax≈b with individual uncertainties not only in b but also in A is an important task. The problem is known in literature as weighted total least squares. In the most general case, correlations between the elements of A,b exist as well. The problem becomes more complicated and can—except for special cases—only be solved numerically. While the formulation of this problem and even its solution is straightforward, its implementation—when the focus is on reliability and computational costs—is not. In this paper, a robust, fast, and universal method for computing the solution of such linear systems as well as their covariance matrix is presented. The results were confirmed by applying the method to several special cases for which an analytical or numerical solution is available. If individual coefficients can be considered to be free of errors, this can be taken into account in a simple way. An implementation of the code in MATLAB is provided.


Introduction
In our laboratory at the National Metrology Institute in Germany, we work in the field of Mueller matrix ellipsometry, an optical method for the investigation of structured surfaces and thin films. While evaluating the raw data (incl. error propagation), we encountered a mathematical problem whose technical solution we think could also be relevant for other experiments.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
It is a classic problem: let A ∈ R m×n with m > n be a coefficient matrix with full column-rank. b ∈ R m×1 and x ∈ R n×1 are column vectors. Typically, the rows of A and b result from m different measurements performed under different configurations or conditions. x is unknown; it will be determined from The solution of this linear system is the central topic of this paper. The solution, its uncertainty, and also the complexity of its calculation strongly depend on the coefficients that are subject to uncertainties, and how these uncertainties are correlated with each other. One searches, therefore, for the bestfit solution x, considering the covariances of the matrix [A, b] and wants to know how these input covariances propagate to W (m·(n+1)×m·(n+1)) = Σ −1 [A,b] . (2) With the vectorization of a matrix performed column-wise the problem can be formulated as an optimization problem: with (5) s.t.
(A + dA) x = b + db (6) ( T denotes the matrix transpose operation; SE, squared error). Or, with the equality constraints eliminated: This is the most general formulation of the problem. Depending on the structure of W, which can be partitioned as (8) many different solutions were found.

Cases with analytical solutions
Analytical solutions were found for . . , m · n 1 , j = m · (n 1 + 1) , . . . , m · n, and n 1 the subset of columns in A, which are assumed to be free of errors. • GTLS (generalized total least squares): W ∼ Here, I is the identity matrix, ⊗ the Kronecker product, and diag (v) generates a diagonal matrix from a vector v. An overview of the methods is given in [1]. A large collection of references can be found in [2].
In the appendix, we provide in brief the formulas for some analytical solutions.

Cases with numerical solutions only
For other cases, however, only numerical solutions are available. These were often developed in the field of geodesy.
Premoli and Rastello have dealt with the element-wise weighted total least squares [3]. This method allows independent weights for the elements of [A, b] but no correlations i.e. W is a diagonal matrix with independent diagonal elements. Later, this method was taken up by further authors and adapted for special applications [4][5][6].
For the case of no correlation between A and b (i.e. W Ab = 0), Schaffrin and Wieser [7] solved the problem for a certain structure of W A : it must be a Kronecker product of a m × m with a n × n matrix. Mahboub [8] succeeded in removing this restriction. Amiri-Simkooei and Jazaeri showed that this weighted total least squares (WTLS) problem is an extension of WLS problem [9].
With the above methods, it is possible to fit a straight line on 2D data. This is what Krystek and Anton did in [10]. Their clou is that they transferred the problem to polar coordinates and reduced it to only one free parameter, which makes their solution very robust. They later modified their method and allowed mutual correlations [11]. Here, W Ab now comes into play and is no longer equal to 0. In addition, in both implementations, they also provided a numerical estimation for the covariance matrix of the solution.
Snow [12] and Zhou et al [13] provided a solution for the case of a full weight matrix W. Later, Zhou and Fang [14] also provided a solution for the mixed case (WLS-WTLS), in which some columns of A are assumed to be free of errors.
Despite the large number of different methods, we must conclude that our problem 1 (and certainly other problems as well) cannot be solved by any of the above methods on account of individual limitations. We have the following requirements: (1) W must be allowed to be a full matrix i.e. correlations between A and b must be allowed. (2) In some of our m measurements, some non-error-free components are not involved. This means it must be possible to exclude individual elements from the WTLS optimization and not only complete columns, like in WLS-WTLS or MTLS.
(3) The method should be fast since many modern problems in metrological practice deal with large matrices, and often a data preprocessing is required during a measurement to control the experimental parameters for the next measurement step 1 . (4) To allow an error-propagation analysis, a covariance matrix for the best fitting solution x must be provided.
The last point especially is of the utmost importance in metrology. In our application, the solution of the problem is an intermediate step; in terms of traceability, we must be able to determine the influence of input uncertainties on the final result and its uncertainty in accordance with the 'Guide to the expression of uncertainty in measurement' [16].

Numerical solution for the general case
In the general case, the weight matrix W has no special structure. But to represent a metric space, it must, of course, be symmetric and positive definite. If it is only positive semi-definite, the problem can be reformulated to only consider the relevant dimensions.
The task is now to solve equation (4). The unknowns are x and dA. For a given x, the necessary condition for an optimal dA is This is a linear system, with m · n equations for the m · n unknowns of da. After exploiting and applying it can be converted to with and da in equation (11) can be solved for a given x in MATLAB with the so-called backslash operator '\' as With da reshaped to dA, one gets SE W from equation (7). SE W then can be used as merit function with a solver of your choice to find the optimal vector x. However, in this form, the optimization is very expensive in terms of computing time and memory requirements. We will address this issue in section 3.

Reduction to requested elements
If individual elements of da and/or db are to be excluded from optimization (i.e. it is assumed to be free of errors), this can be done very easily: one has to just eliminate the regarding rows and columns in the matrices and vectors in equation (11).
For this purpose, let dAi be given as a user-defined Boolean array, indicating the matrix elements of dA to be optimized. dbi is accordingly defined. Furthermore, elements of dAi and dbi shall be set to false by default, if the according diagonal element of the weight matrix W equals 0. Then, equation (11) has to be solved for the subset da dai of da with (15) ('red' = reduced; logical indexing is used here: [T, F, T,…]→[1, 3,…]). Of course, db has to be again calculated from equation (6) taking care that false elements of dAi lead to zeros at the respective positions in dA. SE W (x) can then be calculated from da dai , db dbi , and W red . Note that this technique is more general than MTLS or WLS-WTLS since it allows one to assume individual elements of [A, b] to be free of errors and not only full columns.

Performance and robustness
In the previous chapter, the mathematical solution of the problem was given. We will now cover some technical and programming aspects to make it fast and reliable. The method was implemented in MATLAB. Nevertheless, most of the following aspects can certainly be transferred to other programming languages.

Recommended optimizer
Since the merit function (equation (7)) gives back a scalar positive value-which is the sum of SEs in the metric space defined by W 2 -any multivariate optimizer should be applicable.
We recommend MATLAB's optimizer lsqnonlin, which is part of the optimization toolbox. Its special feature: it expects the merit function to provide an error vector F instead of the sum of squared errors F 2 i . Here, such a vector can easily be defined: (16) with R T W R W = W is the Cholesky decomposition of W. It is obvious that F contains more details on the problem than SE W = F 2 i . This improves convergence and therefore performance. lsqnonlin uses the trust-regionreflective algorithm by default and alternatively the Levenberg-Marquardt. It is a gradient-based local solver. If your problem is difficult (i.e. many local minima), you may need to use a global solver (like differential evolution or particle swarm). In that case, there is nothing left but to use

Initial values for optimization
For good conditioned problems, it is easy and cheap to provide a good guess as the initial value for the optimization from the LS solution (equation (40)) or the WLS solution (equation (41)). If the problem is close to the TLS problem, one should, of course, start with its solution as the initial guess (equation (43)). However, GTLS is the analytical solution that provides the widest possible consideration of a general weight matrix W and should, therefore, be preferred in general.
To do so, first, the necessary matrices P C and P R have to be derived as the solution from the minimization problem: (∥·∥ F denotes the Frobenius norm). This is known as the nearest Kronecker product problem and has an analytical solution [17].
Let R W −1 be a certain rearrangement 3 of W −1 depending on the dimensions of P C and P R , i.e. (n + 1) × (n + 1) and m × m. From its singular value decomposition P C and P R can be found from vec (P C ) = σ 1 U :,1 and vec (P R ) = V :,1 respectively. It is important to note that only the product, P C ⊗ P R can be uniquely determined. P C and P R are unique except for constant scalar factors s and 1/s respectively. Furthermore, since W is symmetric and positive definite, P C and P R are symmetric, and either both positive or both negative definite. In the latter case, they can be made positive definite with s = −1. The absolute value of s has no impact on the GTLS solution, which then gives the best possible analytical guess as initial values for the subsequent optimization. (11) The essential part of the solution given in the previous chapter is equation (11). The terms of this formal expression can be rewritten in a much more memory-(and, hence, time-) efficient way as (19) and

Robust solution of equation (11)
Since W is symmetric, and positive definite, the same is true for the coefficient matrix V T WV . Because of that, equation (11) can be solved efficiently and robust. Let R with be the Cholesky decomposition of the coefficient matrix. Because of R's triangular shape, equation (11) can be solved with a forward-, followed by a backward-substitution. Again, with MATLAB's backslash operator, one can write these two steps as: This solution looks more complicated than the one given in equation (14). But, actually, it is not; there are two reasons to solve da with equation (22) instead of equation (14).
(a) Applying the backslash-operator causes MATLAB to check the coefficient matrix for certain properties. If it is of triangular shape, then the substitution technique is directly applied, and all other tests are skipped. On the other hand, checking whether the Cholesky solver is applicable to equation (14) (which we already know) costs three more tests on the input matrix, which can be saved here. (b) We explicitly will need the matrix R in the further process.

Calculation of the Jacobian matrix
As stated earlier, MATLAB's lsqnonlin is gradient-based. This means that after each calculation of a candidate solution, the Jacobian matrix is calculated to estimate the candidate of the subsequent iteration. Since the accuracy requirements on the Jacobian are relaxed, it is sufficient to estimate it from forward finite differences, which can be obtained after n explicit calls of the vectorial merit function (equation (16)) i.e.: The vector ∆x i has zero elements, barring the ith element, which equals ∆x i . Of course, its calculation can be very time-consuming even after implementing the improvements suggested in the last two sections. The most costly one is the calculation of R, the Cholesky decomposition. It requires ∼O (m · n) 3 operations. Hence, with the forward finite-difference estimation of the Jacobian, n additional Cholesky decompositions would be necessary (∼ O (m · n) 3 · n ).
But this can be avoided. MATLAB's lsqnonlin can also be supported with an externally provided Jacobian. We will now calculate it with only ∼O (m · n) 2 · n operations.
The finite difference method will also be used now. At first, one may try to find a cheap solution to update R when an element of x updates slightly from x i → x i + ∆x i . Only one method is known for this: the rank-one update. Unfortunately, the conditions for applicability are not fulfilled for this here because, with equation (19), it can be shown that an x-update causes a rank-m update for n = 1 and a rank-2m update for n > 1.
Here is the alternative solution: first, let us define shortcuts for the expressions in equation (11) (or equations (19) and (20)) Suppose we had updated x → x + ∆x i . Formally, da(x + ∆x i ) can then be calculated as With the following is surely true: On taking a look at equation (19), it is revealed that dY (∆x i ) is-related to the spectral radius-surely small compared to Y (x) (keep in mind that x = 0 is never a solution of your problem). This allows a first-order Taylor approximation for the matrix inverse: With the shortcut one gets And finally Note: it was not necessary to explicitly calculate the inverse of Y. One can now reshape ∆da (∆x i ) to ∆dA (∆x i ) and calculate The Jacobian can then be estimated as Remark: for this forward difference approximation of the Jacobian, a step size of is used with the constant ε = 2 −52 , which is the distance from 1.0 to the next larger double-precision number, and sign (x) = 1 for x ⩾ 0 and sign (x) = −1 for x < 0. This definition of ∆x i is based on the default definition for this approximation in MATLAB.
In sum, these five steps carried out in this chapter drastically accelerate the optimization. If one has a supported graphics processing unit (GPU) and MATLAB's parallel-computing toolbox, one can even increase the speed further, since many matrix operations used here can benefit from the parallel organization of a GPU. For this purpose, only the input matrices must be defined as gpuArrays. A parallel distribution to several cores of the CPU-on the other hand-cannot be sensibly implemented, as it has no positive effect on account of overhead.
We finish this chapter with some numbers on the performance: table 1 shows a small benchmark for two different machines with or without GPU support and with or without the Jacobian provided. In our application, the problem with A (140×15) is typically good conditioned, so a few ten iterations will suffice. A complete optimization, including the calculation of the covariance matrix (see next chapter) of the solution, is, therefore, finished after 2.5 s at the latest. Since we typically measure data which lead to several thousand such problems, it makes a big difference in terms of applicability and feasibility whether their solution takes a week or only a few hours.

Covariance matrix of the solution
As a metrologist, the writer is not only interested in the final result x but also in its statistical uncertainty. Therefore, the calculation of the covariance matrix Σ x was also implemented. One oft-used estimation of the covariance matrix calculated from the Jacobian at the optimum of a least-squares-merit function is Or, if the weight matrix is known only up to an unknown factor: Since all values are available at the end of the optimization, it can be directly calculated. But while this estimation might be fine for low-dimensional problems it gets worse with increasing n.
A more reliable approximation of the covariance matrix is based on the Hessian matrix H: Or, again, if the weights are only known relative to each other: The technique applied to determine the Jacobian can also be used to calculate the Hessian. Its elements are now calculated with the central finite differences approximation. This is more costly than the forward finite differences approximation but more accurate: Here, a step size of ∆x i = sign ( Since the Hessian is symmetric, it is sufficient to calculate the n · (n + 1) /2 independent elements only. But as can been seen from equation (39), it requires four explicit calculations of SE W for each element. Hence, the effort is already considerable. Therefore, the calculation of the covariance matrix was implemented as optional output. The same is true for the condition number of V T WV. It might give you some helpful information. But, of course, it requires an additional singular value decomposition.
Before starting an optimization, one might also be interested in the condition number or the efficiency [18] of A.

Applications
The method was implemented in MATLAB as a function, which we named cwtls (correlated weights total least squares). It outputs x as solution of equation (4) and-among others-its covariance Σ ′ x , according to equation (38) (to be compatible with other MATLAB functions). From additional output parameters, the covariance can be calculated, according to equation (37). Table 2 shows some application examples for cwtls and a comparison with-mostly analytical-alternatives. If the optimization has converged to the global minimum (which is the case for the most real-world problems, i.e. not artificially ill-conditioned, and with reasonable starting values), then, in each case, the cwtls solution agrees with the alternative solution-of course-up to the tolerance set as the stopping criterion of the optimization. Also, the covariance matrices provided by cwtls agreed with their alternative calculations. In addition, they were validated with those resulting from the numerically determined Hessian matrix at the optimum. For this purpose, the code from [19] was used.
Since, in the literature, the definitions of the covariance matrices are sometimes somewhat vague or imprecise, the corresponding formulas used here are also given.
Examples of numerical applications are given in the supplementary material (available online at stacks.iop.org/MST/ 33/015017/mmedia) to document agreement of cwtls results with the alternative solutions shown in table 2. Table 2 can also be used as a guide, demonstrating how to prepare the inputs to solve a problem with cwtls.
A (perhaps trivial) remark for the case that an analytical solution exists: if one has absolute knowledge of the weight matrix (i.e. there is no need to rescale Σ x with mse), then one can calculate Σ x -according to the formulas in the tablewithout (or before) knowing x. This should be kept in mind when designing an experiment that focuses on low uncertainties for x.
Please note: cwtls is not intended to replace existing and specialized solutions. Rather, it can be used when specialized solutions do not exist.    The TLS solution is given in equation (43). An alternative formulation is: x

Summary and conclusion
A method to solve overdetermined linear systems with arbitrary correlations of the coefficients was introduced. It was optimized in mathematical and technical aspects in terms of computational speed and robustness. It is flexible in application since it allows one to assume single coefficients to be free of error. The covariance of the solutions is also calculated. The solution itself as well as its covariance were validated by comparison with alternative solutions for special cases. The comparison covered all possible correlation scenarios. From that, it can be concluded that the method works reliably and quite fast for all sufficiently well-conditioned realworld problems.
The MATLAB code, with examples, can be found in the supplementary material. It is intended to make it available under the BSD license [27] at the MATLAB central file exchange.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix. Special cases with analytical solutions
Known analytical solutions for some special cases are given below in their most compact form. Please note: again, only the formulas for the most important cases of well-defined problems are given. That is, it is assumed that all specified inverse matrices exist.

Least squares (LS) solution
For W b ∼ I, W A = 0, and W Ab = 0, so, with only deviations in b but not in A and no correlations allowed, one simply gets the best fitting least squares result as the solution of the Gaussian normal equation as with A + = A T A −1 A T , the Moore-Penrose pseudo inverse of

Weighted least squares (WLS) solution
In the more general case for W b ∼ Σ −1 b with Σ b , the covariance matrix for b, but still W A = 0, and W Ab = 0 the solution is given by

Total least squares (TLS) solution
With W ∼ I (i.e. uniform weights in A and b and no correlations), a total least squares problem is defined. Its solution can be calculated from the singular value decomposition of [A, b]: as Generalized total least squares (GTLS) solution GTLS can deal with a certain structure of deviations in A and b and row-and column-wise correlations [5,25]. The weight matrix is given as Hence, GTLS assumes that all rows and/or all columns have equal covariance matrices. With the Cholesky factors and their inverse C T C C C = P C C T R C R = P R (45) the GTLS problem can be transformed to an ordinary TLS problem. From its TLS-solution, the GTLS-solution can be determined as Here, W C was partitioned as