Total Least-Squares Collocation: An Optimal Estimation Technique for the EIV-Model with Prior Information

Abstract

In regression analysis, oftentimes a linear (or linearized) Gauss-Markov Model (GMM) is used to describe the relationship between certain unknown parameters and measurements taken to learn about them. As soon as there are more than enough data collected to determine a unique solution for the parameters, an estimation technique needs to be applied such as ‘Least-Squares adjustment’, for instance, which turns out to be optimal under a wide range of criteria. In this context, the matrix connecting the parameters with the observations is considered fully known, and the parameter vector is considered fully unknown. This, however, is not always the reality. Therefore, two modifications of the GMM have been considered, in particular. First, ‘stochastic prior information’ (p. i.) was added on the parameters, thereby creating the – still linear – Random Effects Model (REM) where the optimal determination of the parameters (random effects) is based on ‘Least Squares collocation’, showing higher precision as long as the p. i. was adequate (Wallace test). Secondly, the coefficient matrix was allowed to contain observed elements, thus leading to the – now nonlinear – Errors-In-Variables (EIV) Model. If not using iterative linearization, the optimal estimates for the parameters would be obtained by ‘Total Least Squares adjustment’ and with generally lower, but perhaps more realistic precision. Here the two concepts are combined, thus leading to the (nonlinear) ’EIV-Model with p. i.’, where an optimal estimation (resp. prediction) technique is developed under the name of ‘Total Least-Squares collocation’. At this stage, however, the covariance matrix of the data matrix – in vector form – is still being assumed to show a Kronecker product structure.

1. Introduction

Over the last 50 years or so, the (linearized) Gauss-Markov Model (GMM) as standard model for the estimation of parameters from collected observation [1,2] has been refined in a number of ways. Two of these will be considered in more detail, namely

• the GMM after strengthening the parameters through the introduction of “stochastic prior information”. The relevant model will be the “Random Effects Model (REM)”, and the resulting estimation technique has become known as “least-squares collocation” [3,4].
• the GMM after weakening the coefficient matrix through the replacement of fixed entries by observed data, resulting in the (nonlinear) Errors-In-Variables (EIV) Model. When nonlinear normal equations are formed and subsequently solved by iteration, the resulting estimation technique has been termed “Total Least-Squares (TLS) estimation” [5,6,7]. The alternative approach, based on iteratively linearizing the EIV-Model, will lead to identical estimates of the parameters [8].

After a compact review and comparison of several key formulas (parameter estimates, residuals, variance component estimate, etc.) within the above three models, a new model will be introduced that allows the strengthening of the parameters and the weakening of the coefficient matrix at the same time. The corresponding estimation technique will be called “TLS collocation” and follows essentially the outline that had first been presented by this author in June 2009 at the Intl. Workshop on Matrices and Statistics in Smolenice Castle (Slovakia); cf. Schaffrin [9].

Since then, further computational progress has been made, e.g. Snow and Schaffrin [10]; but several open questions remain that need to be addressed elsewhere. These include issues related to a rigorous error propagation of the “TLS collocation” results; progress may be achieved here along similar lines as in Snow and Schaffrin [11], but is beyond the scope of the present paper.

For a general overview of the participating models, the following diagram (Figure 1) may be helpful.

The most informative model defines the top position, and the least informative model is at the bottom. The new model can be formed at the intermediate level (like the GMM), but belongs to the “nonlinear world” where nonlinear normal equations need to be formed and subsequently solved by iteration.

2. A Compact Review of the “Linear World”

2.1. The (linearized) Gauss-Markov Model (GMM)

Let the Gauss-Markov Model be defined by

$$y=\underset{n\times m}{A}\xi +e ,\hspace{1em}q:=\mathrm{rk} A\le \min \{m,n\} ,\hspace{1em}e\sim (0,{\sigma }_0^2{I}_n) ,$$

(1)

possibly after linearization and homogeneization; cf. Koch [2], or Grafarend and Schaffrin [12] among many others.

Here, $y$ denotes the $n\times 1$ vector of (incremental) observations,

$A$

the $n\times m$ matrix of coefficients (given),

$\xi$

the $m\times 1$ vector of (incremental) parameters (unknown),

$e$

the $n\times 1$ vector of random errors (unknown) with expectation $E\left\{e\right\}=0$ while the dispersion matrix $D\left\{e\right\}={\sigma }_0^2{I}_n$ is split into the (unknown) factor

${\sigma }_0^2$

as variance component (unit-free) and

$I_n$

as (homogeneized) $n\times n$ symmetric and positive-definite “cofactor matrix” whose inverse is better known as “weight matrix” $P$; here $P=I_n$ for the sake of simplicity.

Now, it is well known that the Least-Squares Solution (LESS) is based on the principle

$$e^{\mathrm{T}}e=\min \mathrm{s}.\mathrm{t}.\hspace{1em}e=y-A \xi$$

(2)

which leads to the “normal equations”

$$N \widehat{\xi }=c \hspace{1em}\mathrm{for}\hspace{1em}[N, c]:=A^{\mathrm{T}}[A, y].$$

(3)

Depending on rk N, the rank of the matrix N, the LESS may turn out uniquely as

$${\widehat{\xi }}_{\mathrm{LESS}}=N^{-1}c \hspace{1em}\mathrm{iff}\hspace{1em}\mathrm{rk} N=\mathrm{rk} A=q=m;$$

(4)

or it may belong to a solution (hyper)space that can be characterized by certain generalized inverses of N, namely

$$\begin{array}{ll}{\widehat{\xi }}_{\mathrm{LESS}}& \in \left\{N^{-}c\right|N N^{-}N=N\}=\\ & =\left\{N_{rs}^{-}c\right|N N_{rs}^{-}N=N , N_{rs}^{-}N N_{rs}^{-}=N_{rs}^{-}={(N_{rs}^{-})}^{\mathrm{T}}\} ;\end{array}$$

(5)

where $N_{rs}^{-}$ denotes an arbitrary reflexive symmetric g-inverse of N (including the “pseudo-inverse” $N^{+}$),

$$\mathrm{iff}\hspace{1em}\mathrm{rk} N=\mathrm{rk} A=q<m.$$

(6)

In the latter case (6), any LESS will be biased with

$$\mathrm{bias} ({\widehat{\xi }}_{\mathrm{LESS}})=E\{{\widehat{\xi }}_{\mathrm{LESS}}-\xi \}\in \{-(I_m-N^{-}N) \xi \}=\{-(I_m-N_{rs}^{-}N) \xi \},$$

(7)

its dispersion matrix will be

$$D \left\{{\widehat{\xi }}_{\mathrm{LESS}}\right\}={\sigma }_0^2{N}^{-}N{(N^{-})}^{\mathrm{T}}\in \left\{{\sigma }_0^2{N}_{rs}^{-}\right|N N_{rs}^{-}N=N, N_{rs}^{-}N N_{rs}^{-}=N_{rs}^{-}={(N_{rs}^{-})}^{\mathrm{T}}\},$$

(8)

and its Mean Squared Error (MSE) matrix, therefore,

$$MSE \left\{{\widehat{\xi }}_{\mathrm{LESS}}\right\}\in \{{\sigma }_0^2[N_{rs}^{-}+(I_m- N_{rs}^{-}N)(\xi {\sigma }_0^{-2}{\xi }^{\mathrm{T}}){(I_m-N_{rs}^{-}N)}^{\mathrm{T}}]\}.$$

(9)

In contrast to the choice of LESS in case (6), the residual vector will be unique; for any ${\widehat{\xi }}_{\mathrm{LESS}}$:

$${\tilde{e}}_{\mathrm{LESS}}=y-A {\widehat{\xi }}_{\mathrm{LESS}} \sim (0, {\sigma }_0^2(I_n-A N_{rs}^{-}{A}^{\mathrm{T}})=D\{y\}-D\{A {\widehat{\xi }}_{\mathrm{LESS}}\}).$$

(10)

Hence, the optimal variance component estimate will also be unique:

$${\widehat{\sigma }}_0^2={(n-q)}^{-1}{\tilde{e}}_{\mathrm{LESS}}^{\mathrm{T}}\cdot {\tilde{e}}_{\mathrm{LESS}}={(n-q)}^{-1}(y^{\mathrm{T}}y-c^{\mathrm{T}}{\widehat{\xi }}_{\mathrm{LESS}}),$$

(11)

$$E\left\{{\widehat{\sigma }}_0^2\right\}={\sigma }_0^2 , D\left\{{\widehat{\sigma }}_0^2\right\}\approx 2{({\sigma }_0^2)}^2/(n-q) \mathrm{under} \text{quasi-normality}.$$

(12)

For the corresponding formulas in the full-rank case (4), the reflexive symmetric g-inverse $N_{rs}^{-}$ simply needs to be replaced by the regular inverse $N^{-1}$, thereby showing that the LESS turns into an unbiased estimate of $\xi$ while its MSE-matrix coincides with its dispersion matrix (accuracy $\stackrel{\wedge}{=}$ precision).

2.2. The Random Effects Model (REM)

Now, additional prior information (p. i.) is introduced to strengthen the parameter vector within the GMM (1). Based on Harville’s [13] notation of a $m\times 1$ vector $\underset{~}{0}$ of so-called “pseudo-observations”, the new $m\times 1$ vector

$$x:=\xi -\underset{~}{0}\hspace{1em}\hspace{1em}\mathrm{of} “\mathrm{random} \mathrm{effects}” \left(\mathrm{unknown}\right) \mathrm{is} \mathrm{formed} \mathrm{with}$$

(13)

$$E\left\{x\right\}={\beta }_0\hspace{1em}\hspace{1em}\mathrm{as} \mathrm{given} m\times 1 \mathrm{vector} \mathrm{of} \mathrm{p}. \mathrm{i}. \mathrm{and}$$

(14)

$$D\left\{x\right\}={\sigma }_0^2{Q}_0\hspace{1em}\hspace{1em}\mathrm{as} n\times n \text{positive-(semi)definite} \mathrm{dispersion} \mathrm{matrix} \mathrm{of} \mathrm{the} \mathrm{p}. \mathrm{i}.$$

(15)

Consequently, the Random Effects Model (REM) can be stated as

$$\underset{~}{y}:=y-A\cdot \underset{~}{0}=\underset{n\times m}{A}x+e ,\hspace{1em}q:=\mathrm{rk} A\le \min \{m,n\} ,\hspace{1em}\underset{m\times 1}{{\beta }_0}:=x+e_0,$$

(16)

$$\left[\begin{array}{c}e\\ e_0\end{array}\right]\sim ( \left[\begin{array}{c}0\\ 0\end{array}\right] , {\sigma }_0^2\left[\begin{array}{cc}{I}_n& 0\\ 0& Q_0\end{array}\right] ) ,\hspace{1em}\underset{m\times m}{Q_0} \mathrm{symmetric} \mathrm{and} \mathrm{nnd} (\text{non-negative} \mathrm{definite}).$$

(17)

If $P_0:=Q_0^{-1}$ exists, the estimation/prediction of x may now be based on the principle

$$e^{\mathrm{T}}e+e_0^{\mathrm{T}}{P}_0 e_0=\min \hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{1em}e=\underset{~}{y}-Ax ,\hspace{1em}{e}_0={\beta }_0-x,$$

(18)

which leads to the “normal equations” for the “Least-Squares Collocation (LSC)” solution

$$(N+P_0)\cdot \tilde{x}=\underset{~}{c} +P_0\cdot {\beta }_0$$

(19)

which is unique even in the rank-deficient case (6).

If $Q_0$ is singular, but $e_0\in \Re (Q_0)$ with probability 1, then there must exist an $m\times 1$ vector ${\nu }^0$ with

$$Q_0 {\nu }^0=-e_0=x-{\beta }_0\hspace{1em}\mathrm{with} \mathrm{probability} 1 (\mathrm{a}.\mathrm{s}.\stackrel{\wedge}{=}\mathrm{almost} \mathrm{surely}).$$

(20)

Thus, the principle (18) may be equivalently replaced by

$$e^{\mathrm{T}}e+{({\nu }^0)}^{\mathrm{T}}{Q}_0 {\nu }^0=\min \hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{1em}e=\underset{~}{y}-Ax \mathrm{and} (20),$$

(21)

which generates the LSC solution uniquely from the “modified normal equations”

$$\boxed{ (I_m+Q_{0}N)\cdot {\tilde{x}}_{\mathrm{LSC}}={\beta }_0+Q_0\underset{~}{c} } \mathrm{for} \underset{~}{c}:=A^{\mathrm{T}}\underset{~}{y}$$

(22)

or, alternatively, via the “update formula”

$$\boxed{ {\tilde{x}}_{\mathrm{LSC}}={\beta }_0+Q_0 {(I_m+N{Q}_0)}^{-1}(\underset{~}{c}-N{\beta }_0) }$$

(23)

which exhibits the “weak (local) unbiasedness” of ${\tilde{x}}_{\mathrm{LSC}}$ via

$$E\left\{{\tilde{x}}_{\mathrm{LSC}}\right\}={\beta }_0+Q_0 {(I_m+N{Q}_0)}^{-1}(E\{\underset{~}{c}\}-N{\beta }_0)={\beta }_0+0=E\left\{x\right\}.$$

(24)

Consequently, the MSE-matrix of ${\tilde{x}}_{\mathrm{LSC}}$ can be obtained from

$$MSE \left\{{\tilde{x}}_{\mathrm{LSC}}\right\}=D \{{\tilde{x}}_{\mathrm{LSC}}-x\}={\sigma }_0^2{Q}_0 {(I_m+N{Q}_0)}^{-1}=$$

(25)

$$={\sigma }_0^2 {(N+P_0)}^{-1}\hspace{1em}\mathrm{if}\hspace{1em}{P}_0:=Q_0^{-1} \mathrm{exists},$$

(26)

whereas the dispersion matrix of ${\tilde{x}}_{\mathrm{LSC}}$ itself is of no relevance here. It holds:

$$D\left\{{\tilde{x}}_{\mathrm{LSC}}\right\}=D \left\{x\right\}-D \{x-{\tilde{x}}_{\mathrm{LSC}}\}={\sigma }_0^2{Q}_0- {\sigma }_0^2{Q}_0{(I_m+N{Q}_0)}^{-1}.$$

(27)

Again, both residual vectors are uniquely determined from

$${({\tilde{e}}_0)}_{\mathrm{LSC}}={\beta }_0-{\tilde{x}}_{\mathrm{LSC}}=-Q_0{\widehat{\nu }}_{\mathrm{LSC}}^0\hspace{1em}\mathrm{for}$$

(28)

$${\widehat{\nu }}_{\mathrm{LSC}}^0:={(I_m+N{Q}_0)}^{-1}(\underset{~}{c} -N{\beta }_0), \mathrm{and}$$

(29)

$${\tilde{e}}_{\mathrm{LSC}}=\underset{~}{y} -A {\tilde{x}}_{\mathrm{LSC}}=[I_m-A Q_0{(I_m+N{Q}_0)}^{-1}{A}^{\mathrm{T}}] (\underset{~}{y} -A{\beta }_0)$$

(30)

with the control formula

$$\boxed{ A^{\mathrm{T}}{\tilde{e}}_{LSC}={\widehat{\nu }}_{\mathrm{LSC}}^0 } .$$

(31)

An optimal estimate of the variance component may now be obtained from

$$\begin{array}{cc}{({\widehat{\sigma }}_0^2)}_{\mathrm{LSC}}\hfill & =n^{-1}\cdot ({\tilde{e}}_{\mathrm{LSC}}^{\mathrm{T}}\cdot {\tilde{e}}_{\mathrm{LSC}}+{({\widehat{\nu }}_{\mathrm{LSC}}^0)}^{\mathrm{T}}{Q}_0 {\widehat{\nu }}_{\mathrm{LSC}}^0)=\hfill \\ \hfill & =n^{-1}\cdot ({\underset{~}{y}}^{\mathrm{T}}\underset{~}{y} - {\underset{~}{c}}^{\mathrm{T}}{\tilde{x}}_{\mathrm{LSC}}-{\beta }_0^{\mathrm{T}}{\widehat{\nu }}_{\mathrm{LSC}}^0) .\hfill \end{array}$$

(32)

3. An Extension into the “Nonlinear World”

3.1. The Errors-In-Variables (EIV) Model

In this scenario, the Gauss-Markov Model (GMM) is further weakened by allowing some or all of the entries in the coefficient matrix A to be observed. So, after introducing a corresponding $n\times m$ matrix of unknown random errors $E_A$, the original GMM (1) turns into the EIV-Model

$$y=\underset{n\times m}{\underbrace{(A-E_A)}}\cdot \xi +e ,\hspace{1em}q:=\mathrm{rk} A=m<n ,\hspace{1em}e\sim (0,{\sigma }_0^2{I}_n) ,$$

(33)

with the vectorized form of $E_A$ being characterized through

$$\underset{nm\times 1}{e_A}:=\mathrm{vec} E_A\sim (0,{\sigma }_0^2(I_m\otimes I_n)={\sigma }_0^2{I}_{mn}) ,\hspace{1em}C\{e,e_A\}=0 \left(\mathrm{assumed}\right).$$

(34)

Here, the vec operation transform a matrix into a vector by stacking all columns underneath each other while $\otimes$ denotes the Kronecker-Zehfuss product of matrices, defined by

$$\underset{p\times q}{G}\otimes \underset{r\times s}{H}={[g_{ij}\cdot H]}_{pr\times qs}\hspace{1em}\mathrm{if}\hspace{1em}G=[g_{ij}].$$

(35)

In particular, the following key formula holds true:

$$\mathrm{vec} (A B C^{\mathrm{T}})=(C\otimes A)\cdot \mathrm{vec} B$$

(36)

for matrices of suitable size. Note that, in (34), the choice $Q_A:=I_{mn}$ is a very special one. In general, $Q_A$ may turn out singular whenever some parts of the matrix A remain unobserved (i.e., nonrandom).

In any case, thanks to the term

$$E_A\cdot \xi =({\xi }^{\mathrm{T}}\otimes I_n) e_A,$$

(37)

the model (33) needs to be treated in the “nonlinear world” even though the vector $\xi$ may contain only incremental parameters. From now on, A is assumed to have full column rank, $\mathrm{rk} A=:q=m$.

Following Schaffrin and Wieser [7], for instance, the “Total-Least Squares Solution (TLSS)” can be based on the principle

$$e^{\mathrm{T}}e+e_A^{\mathrm{T}}{e}_A=\min \hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{1em}(33\text{-}34),$$

(38)

and the Lagrange function

$$\begin{array}{ll}\Phi (e,e_A,\xi ,\lambda ):& =e^{\mathrm{T}}e+e_A^{\mathrm{T}}{e}_A+2{\lambda }^{\mathrm{T}}(y-A \xi -e+E_A\xi )=\\ & =e^{\mathrm{T}}e+e_A^{\mathrm{T}}{e}_A+2{\lambda }^{\mathrm{T}}[y-A \xi -e+({\xi }^{\mathrm{T}}\otimes I_n) e_A]\end{array}$$

(39)

which needs to be made stationary. The necessary Euler-Lagrange conditions then read:

$$\frac{1}{2}\frac{\partial \Phi }{\partial e}=\tilde{e}-\widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}\widehat{\lambda }=\tilde{e}$$

(40)

$$\frac{1}{2}\frac{\partial \Phi }{\partial e_A}={\tilde{e}}_A+({\widehat{\xi }}_{\mathrm{TLS}}\otimes I_n) \widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}{\tilde{E}}_A=-\widehat{\lambda } {\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}}$$

(41)

$$\frac{1}{2}\frac{\partial \Phi }{\partial \xi }=-{(A-{\tilde{E}}_A)}^{\mathrm{T}}\widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}{A}^{\mathrm{T}}\widehat{\lambda }={\tilde{E}}_A^{\mathrm{T}}\widehat{\lambda }=-{\widehat{\xi }}_{\mathrm{TLS}}\cdot ({\widehat{\lambda }}^{\mathrm{T}}\widehat{\lambda })$$

(42)

$$\frac{1}{2}\frac{\partial \Phi }{\partial \lambda }=y-A {\widehat{\xi }}_{\mathrm{TLS}}-\tilde{e}+{\tilde{E}}_A {\widehat{\xi }}_{\mathrm{TLS}}\doteq 0\hspace{1em}\Rightarrow \hspace{1em}y-A {\widehat{\xi }}_{\mathrm{TLS}}=\widehat{\lambda } (1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})\hspace{1em}$$

(43)

$$\Rightarrow \hspace{1em}c-N{\widehat{\xi }}_{\mathrm{TLS}}:=A^{\mathrm{T}}(y-A {\widehat{\xi }}_{\mathrm{TLS}})=A^{\mathrm{T}}\widehat{\lambda } (1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})=-{\widehat{\xi }}_{\mathrm{TLS}}\cdot {\widehat{\nu }}_{\mathrm{TLS}}$$

(44)

for

$${\widehat{\nu }}_{\mathrm{TLS}}:=({\widehat{\lambda }}^{\mathrm{T}}\widehat{\lambda })\cdot (1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})={\widehat{\lambda }}^{\mathrm{T}}(y-A {\widehat{\xi }}_{\mathrm{TLS}})=$$

(45)

$$\hspace{1em}={(1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})}^{-1}\cdot {(y-A {\widehat{\xi }}_{\mathrm{TLS}})}^{\mathrm{T}}(y-A {\widehat{\xi }}_{\mathrm{TLS}})=$$

(46)

$$={(1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})}^{-1} [y^{\mathrm{T}}(y-A {\widehat{\xi }}_{\mathrm{TLS}})-{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} (c-N {\widehat{\xi }}_{\mathrm{TLS}})] \ge 0$$

(47)

$$\Rightarrow (1+\underset{}{\underbracket{{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})\cdot {\widehat{\nu }}_{\mathrm{TLS}}}}=y^{\mathrm{T}}y-c^{\mathrm{T}}{\widehat{\xi }}_{\mathrm{TLS}}+\underset{}{\underbracket{({\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})\cdot {\widehat{\nu }}_{\mathrm{TLS}}}}$$

$$\Rightarrow \boxed{ {\widehat{\nu }}_{\mathrm{TLS}}=y^{\mathrm{T}}y-c^{\mathrm{T}}{\widehat{\xi }}_{\mathrm{TLS}} }$$

(48)

which needs to be solved in connection with the “modified normal equations” from (44), namely

$$\boxed{ (N-{\widehat{\nu }}_{\mathrm{TLS}}{I}_m) {\widehat{\xi }}_{\mathrm{TLS}}=c } .$$

(49)

Due to the nonlinear nature of ${\widehat{\xi }}_{\mathrm{TLS}}$, it is not so easy to determine if it is an unbiased estimate, or how its MSE-matrix may exactly look like. First attempts of a rigorous error propagation have recently been undertaken by Amiri-Simkooei et al. [14] and by Schaffrin and Snow [15], but are beyond the scope of this paper.

Instead, both the optimal residual vector ${\tilde{e}}_{\mathrm{TLS}}$ and the optimal residual matrix ${({\tilde{E}}_A)}_{\mathrm{TLS}}$ are readily available through (40) and (43) as

$${\tilde{e}}_{\mathrm{TLS}}=\widehat{\lambda }=(y-A {\widehat{\xi }}_{\mathrm{TLS}})\cdot {(1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})}^{-1},$$

(50)

and through (41) as

$${({\tilde{E}}_A)}_{\mathrm{TLS}}=-{\tilde{e}}_{\mathrm{TLS}}\cdot {\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}}=-(y-A {\widehat{\xi }}_{\mathrm{TLS}})\cdot {(1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})}^{-1}\cdot {\widehat{\xi }}_{\mathrm{TLS}}.$$

(51)

As optimal variance component estimate, it is now proposed to use the formula

$${({\widehat{\sigma }}_0^2)}_{\mathrm{TLS}}={(n-m)}^{-1}\cdot [ {\tilde{e}}_{\mathrm{TLS}}^{\mathrm{T}} {\tilde{e}}_{\mathrm{TLS}}+{({\tilde{e}}_A)}_{\mathrm{TLS}}^{\mathrm{T}} {({\tilde{e}}_A)}_{\mathrm{TLS}}]=$$

$$\hspace{1em} ={(n-m)}^{-1}\cdot {\widehat{\lambda }}^{\mathrm{T}}\widehat{\lambda }(1+{\widehat{\xi }}_{\mathrm{TLS}}^{\mathrm{T}} {\widehat{\xi }}_{\mathrm{TLS}})={\widehat{\nu }}_{\mathrm{TLS}}/(n-m),$$

(52)

in analogy to the previous estimates (11) and (32).

3.2. A New Model: The EIV-Model with Random Effects (EIV-REM)

In the following, the above EIV-Model (33-34) is strengthened by introducing stochastic prior information (p. i.) on the parameters which thereby change their character and become “random effects” as in (13-15). The EIV-REM can, therefore, be stated as

$$\underset{~}{y}=\underset{n\times m}{\underbrace{(A-E_A)}}\cdot x+e ,\hspace{1em}q:=\mathrm{rk} A\le \min \{m,n\} ,\hspace{1em}{\beta }_0=x+e_0 (\mathrm{given}),$$

(53)

with

$$\left[\begin{array}{c}e\\ e_A=\mathrm{vec} E_A\\ e_0\end{array}\right] \sim ( \left[\begin{array}{c}0\\ 0\\ 0\end{array}\right] , {\sigma }_0^2\left[\begin{array}{ccc}{I}_n& 0& 0\\ 0& I_{mn}& 0\\ 0& 0& Q_0\end{array}\right] ) , Q_0 \mathrm{symmetric} \mathrm{and} \mathrm{nnd}.$$

(54)

The first set of formulas will be derived by assuming that the weight matrix $P_0:=Q_0^{-1}$ exists uniquely for the p. i. Then, the “TLS collocation (TLSC)” may be based on the principle

$$e^{\mathrm{T}}e+e_A^{\mathrm{T}}{e}_A+e_0^{\mathrm{T}}{P}_0 e_0=\min \hspace{1em}\mathrm{s}.\mathrm{t}.\hspace{1em}(53\text{-}54),$$

(55)

resp. on the Lagrange function

$$\Phi (e,e_A,e_0^{},\lambda ):=e^{\mathrm{T}}e+e_A^{\mathrm{T}}{e}_A+e_0^{\mathrm{T}}{P}_0 e_0+2{\lambda }^{\mathrm{T}}[(\underset{~}{y} -A {\beta }_0-e+({\beta }_0^{\mathrm{T}}\otimes I_n) e_A+A e_0\underset{=-(e_0^{\mathrm{T}}\otimes I_n)\cdot e_A}{\underbrace{-E_A e_0}}]$$

(56)

which needs to be made stationary. The necessary Euler-Lagrange conditions then read:

$$\frac{1}{2}\frac{\partial \Phi }{\partial e}=\tilde{e}-\widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}\widehat{\lambda }=\tilde{e}$$

(57)

$$\frac{1}{2}\frac{\partial \Phi }{\partial e_A}={\tilde{e}}_A+[({\beta }_0-{\tilde{e}}_0)\otimes I_n] \widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}{\tilde{E}}_A=-\widehat{\lambda }\cdot {({\beta }_0-{\tilde{e}}_0)}^{\mathrm{T}}=:-\widehat{\lambda }\cdot {\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}}$$

(58)

$$\frac{1}{2}\frac{\partial \Phi }{\partial e_0}=P_0 {\tilde{e}}_0+{(A-{\tilde{E}}_A)}^{\mathrm{T}}\widehat{\lambda }\doteq 0\hspace{1em}\Rightarrow \hspace{1em}{A}^{\mathrm{T}}\widehat{\lambda }={\tilde{E}}_A^{\mathrm{T}}\widehat{\lambda }-P_0 {\tilde{e}}_0\hspace{1em}\Rightarrow \hspace{1em}$$

(59)

$$\Rightarrow \hspace{1em}{A}^{\mathrm{T}}\widehat{\lambda }=-{\tilde{x}}_{\mathrm{TLSC}}\cdot ({\widehat{\lambda }}^{\mathrm{T}}\widehat{\lambda })+{\widehat{\nu }}_{\mathrm{TLSC}}^0\hspace{1em}\mathrm{for}\hspace{1em}{\widehat{\nu }}_{\mathrm{TLSC}}^0:=-P_0 {\tilde{e}}_0=P_0({\beta }_0-{\tilde{x}}_{\mathrm{TLSC}})$$

(60)

$$\frac{1}{2}\frac{\partial \Phi }{\partial \lambda }=\underset{~}{y}-A ({\beta }_0-{\tilde{e}}_0)-\tilde{e}+{\tilde{E}}_A ({\beta }_0-{\tilde{e}}_0)\doteq 0\hspace{1em}\Rightarrow \hspace{1em}\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}}=\widehat{\lambda } (1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})\hspace{1em}\Rightarrow \hspace{1em}$$

(61)

$$\Rightarrow \hspace{1em}\widehat{\lambda }=(\underset{~}{y}-A {\tilde{x}}_{\mathrm{TLSC}})\cdot {(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-1}.$$

(62)

Combining (60) with (62) results in

$$(\underset{~}{c}-N {\tilde{x}}_{\mathrm{TLSC}})\cdot {(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-1}=A^{\mathrm{T}}\widehat{\lambda }=-{\tilde{x}}_{\mathrm{TLSC}}{(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-2}{(\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}})}^{\mathrm{T}}(\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}})+{\widehat{\nu }}_{\mathrm{TLSC}}^0,$$

(63)

and finally in

$$\boxed{ (N+(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})P_0-{\widehat{\nu }}_{\mathrm{TLSC}}{I}_m) {\tilde{x}}_{\mathrm{TLSC}}=\underset{~}{c} +P_0 {\beta }_0\cdot (1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}}) }$$

(64)

where

$$\begin{gathered} \begin{array}{}\boxed{ {\widehat{\nu }}_{\mathrm{TLSC}}:={(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-1}{(\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}})}^{\mathrm{T}}(\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}}) , \mathrm{and} \\ \\ {\widehat{\nu }}_{\mathrm{TLSC}}^0:=-P_0 ({\beta }_0-{\tilde{x}}_{\mathrm{TLSC}})=-P_0 {\tilde{e}}_0 , \mathrm{provided} \mathrm{that} P_0 \mathrm{exists} .}\end{array} \end{gathered}$$

(65)

In the more general case of a singular matrix $Q_0$, an approach similar to (20) can be followed, leading to the equation system

$$\boxed{ [(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})\cdot I_m+Q_{0}N-{\widehat{\nu }}_{\mathrm{TLSC}}\cdot Q_0] {\tilde{x}}_{\mathrm{TLSC}}={\beta }_0\cdot (1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})+Q_0\underset{~}{c} }$$

(66)

that needs to be solved in connection with (65). Obviously,

$$\begin{array}{}\boxed{ {\tilde{x}}_{\mathrm{TLSC}}\hspace{1em}\mapsto {\beta }_0\mathrm{if}\hspace{1em}{Q}_0\hspace{1em}\mapsto \hspace{1em}0 ,\hspace{1em}\mathrm{and} \\ {\tilde{x}}_{\mathrm{TLSC}}\hspace{1em}\mapsto {\tilde{x}}_{\mathrm{TLS}}\hspace{1em} \mathrm{if}\hspace{1em}{P}_0\hspace{1em} \mapsto \hspace{1em}0 .}\end{array}$$

(67)

Again, it is still unclear if ${\tilde{x}}_{\mathrm{TLSC}}$ represents an unbiased prediction of the vector $x$ of random effects. Also, very little (if anything) is known about the corresponding MSE-matrix of ${\tilde{x}}_{\mathrm{TLSC}}$. The answers to these open problems will be left for a future contribution. It is, however, possible to find the respective residual vectors/matrices represented as follows:

$${\tilde{e}}_{\mathrm{TLSC}}=\widehat{\lambda }=(\underset{~}{y}-A {\tilde{x}}_{\mathrm{TLSC}})\cdot {(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-1},$$

(68)

$${({\tilde{E}}_A)}_{\mathrm{TLSC}}=-\widehat{\lambda }\cdot {\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}}=-(\underset{~}{y} -A {\tilde{x}}_{\mathrm{TLSC}})\cdot {(1+{\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}} {\tilde{x}}_{\mathrm{TLSC}})}^{-1}\cdot {\tilde{x}}_{\mathrm{TLSC}}^{\mathrm{T}},$$

(69)

$${({\tilde{e}}_0)}_{\mathrm{TLSC}}=-Q_0\cdot {\widehat{\nu }}_{\mathrm{TLSC}}^0={\beta }_0- {\tilde{x}}_{\mathrm{TLSC}},$$

(70)

while a suitable formula for the variance component is suggested as

$${({\widehat{\sigma }}_0^2)}_{\mathrm{TLSC}}=n^{-1}\cdot [ {\widehat{\nu }}_{\mathrm{TLSC}}+{({\widehat{\nu }}_{\mathrm{TLSC}}^0)}^{\mathrm{T}}{Q}_0 ({\widehat{\nu }}_{\mathrm{TLSC}}^0)].$$

(71)

4. Conclusions and Outlook

Key formulas have been developed successfully to optimally determine the parameters and residuals within the new ‘EIV-Model with p. i.’ (or EIV-REM) which turns out to be more general than the other three models considered here (GMM, REM, EIV-Model). In particular, it is quite obvious that

• EIV-REM becomes the REM if $D\left\{e_A\right\}:=0$,
• EIV-REM becomes the EIV-Model if $P_0:=0$,
• EIV-REM becomes the GMM if both $P_0:=0$ and $D\left\{e_A\right\}:=0$.

Hence the new EIV-REM can indeed serve as a universal representative of the whole class of models presented here.

Therefore, in a follow-up paper, it is planned to also cover more general dispersion matrices for $e$ and $e_A$ in (54), similarly to the work by Schaffrin et al. [16] for the EIV-Model with singular dispersion matrices for $e_A$.