A Low Rank Gaussian Process Prediction Model for Very Large Datasets

TL;DR

This paper introduces an efficient iterative Fixed Rank Kriging algorithm that guarantees positive definiteness of the covariance matrix, enabling scalable spatial prediction for very large datasets, demonstrated on ocean chlorophyll data.

Contribution

It proposes a novel iterative Fixed Rank Kriging method ensuring positive definiteness, improving scalability and reliability for large spatial datasets.

Findings

Algorithm converges under mild conditions.

Predicts missing ocean chlorophyll data effectively.

Outperforms traditional spatial prediction methods.

Abstract

Spatial prediction requires expensive computation to invert the spatial covariance matrix it depends on and also has considerable storage needs. This work concentrates on computationally efficient algorithms for prediction using very large datasets. A recent prediction model for spatial data known as Fixed Rank Kriging is much faster than the kriging and can be easily implemented with less assumptions about the process. However, Fixed Rank Kriging requires the estimation of a matrix which must be positive definite and the original estimation procedure cannot guarantee this property. We present a result that shows when a matrix subtraction of a given form will give a positive definite matrix. Motivated by this result, we present an iterative Fixed Rank Kriging algorithm that ensures positive definiteness of the matrix required for prediction and show that under mild conditions the algorithm numerically converges. The modified Fixed Rank Kriging procedure is implemented to predict missing chlorophyll observations for very large regions of ocean color. Predictions are compared to those made by other well known methods of spatial prediction.

Tables

Table 1. TABLE I: A side to side comparison of some properties of the two ocean color measuring instruments SeaWiFS and MODIS. Adapted from [ 19 ] and http://www.ioccg.org/sensors/current.html.

Sensor	SeaWiFS	MODIS
Satellite	Orbview-2	Aqua
Data available since	09/04/1997	07/04/2002
Time equator is crossed	12:20pm	1:30pm
Spectral bands	8	36
Swath width	2806 km	2330 km
Spectral coverage (nm)	402-885	405-14,385
Resolution	1100 m	1000 m
Tilt	$-{20}^o,0^o,{20}^o$	None
Orbit	Descending	Ascending

Table 2. TABLE II: Comparison of the mean squared prediction error obtained from OLS, AM and FRK models for the region in the North Pacific when 15, 25 and 50% of the data are removed and used as test data. Each cell contains the mean M ​ S ​ P ​ E ^ ^ 𝑀 𝑆 𝑃 𝐸 \widehat{MSPE} for all 50 simulated test data and the corresponding standard deviation in parenthesis.

$𝑴𝒐𝒅𝒆𝒍$	$\widehat{𝑴𝑺𝑷𝑬}$ (15%)	$\widehat{𝑴𝑺𝑷𝑬}$ (25%)	$\widehat{𝑴𝑺𝑷𝑬}$ (50%)	CPU time
OLS	0.0515 (6.50e-4)	0.0517 (5.75e-4)	0.0516 (3.63e-4)	0.27
AM	0.0169 (3.54e-4)	0.0169 (2.71e-4)	0.0169 (1.93e-4)	140.84
FRK	0.0099 (2.73e-4)	0.0100 (1.88e-4)	0.0100 (8.73e-5)	167.68

Table 3. TABLE III: Comparison of the mean squared prediction error obtained from OLS, AM and FRK models for the region in the South Pacific when 15, 25 and 50% of the data are removed and used as test data. Each cell contains the mean M ​ S ​ P ​ E ^ ^ 𝑀 𝑆 𝑃 𝐸 \widehat{MSPE} for all 50 simulated test data and the corresponding standard deviation in parenthesis.

$𝑴𝒐𝒅𝒆𝒍$	$\widehat{𝑴𝑺𝑷𝑬}$ (15%)	$\widehat{𝑴𝑺𝑷𝑬}$ (25%)	$\widehat{𝑴𝑺𝑷𝑬}$ (50%)
OLS	0.0598 (8.00e-4)	0.0597 (5.45e-4)	0.0597 (4.21e-4)
AM	0.0188 (7.33e-4)	0.0188 (5.16e-4)	0.0188 (3.04e-4)
FRK	0.0159 (6.83e-4)	0.0159 (4.73e-4)	0.0159 (2.63e-4)

Equations

$$\displaystyle Y(\mbox{\boldmath$s$})$$

$$\displaystyle C(W(\mbox{\boldmath$s$}_{i}),W(\mbox{\boldmath$s$}_{j}))=\mbox{\boldmath$Z$}(\mbox{\boldmath$s$}_{i})^{{}^{\prime}}K\mbox{\boldmath$Z$}(\mbox{\boldmath$s$}_{j})$$

$$\displaystyle Var_{K}(\mbox{\boldmath$Y$})\equiv\Sigma_{K}=ZKZ^{\prime}+\sigma^{2}V.$$

$${\Sigma }_K^{-1}=({\sigma }^{2}V{)}^{-1}-$$

$$\displaystyle\hat{H}(\mbox{\boldmath$s$}_{o})=\mbox{\boldmath$x$}(\mbox{\boldmath$s$}_{o})^{{}^{\prime}}\mbox{\boldmath$\hat{\beta}$}+\mbox{\boldmath$Z$}(\mbox{\boldmath$s$}_{o})^{\prime}KZ^{\prime}\Sigma^{-1}_{K}(\mbox{\boldmath$Y$}-X\mbox{\boldmath$\hat{\beta}$})$$

$$\displaystyle\sigma^{2}_{K}(\mbox{\boldmath$s$}_{o})$$

$$\displaystyle w_{mi}=\left\{\begin{array}[]{ll}1&\mathrm{if}\hskip 2.84544pt\mbox{\boldmath$s$}_{i}\in N(\mbox{\boldmath$u$}_{m}),\ \mathrm{and}\ Y(\mbox{\boldmath$s$}_{i})\ \mathrm{is}\ \mathrm{not}\ \mathrm{missing},\vspace{0.5cm}\\ 0&\hskip 2.84544pt\rm{otherwise}\end{array}\right.$$

$$\displaystyle\bar{D}(\mbox{\boldmath$u$}_{m})=\frac{\mbox{\boldmath$w$}_{m}^{{}^{\prime}}\mbox{\boldmath$D$}}{\mbox{\boldmath$w$}_{m}^{{}^{\prime}}\mbox{\boldmath$1$}_{n}}\hskip 14.22636pt\mathrm{for}\hskip 14.22636ptm=1,...,M$$

$$\displaystyle\widehat{\Sigma}_{M}(\mbox{\boldmath$u$}_{m},\mbox{\boldmath$u$}_{k})=\left\{\begin{array}[]{ll}V_{D}(\mbox{\boldmath$u$}_{k})&\mathrm{if}\hskip 2.84544ptm=k,\vspace{0.5cm}\\ C_{D}(\mbox{\boldmath$u$}_{m},\mbox{\boldmath$u$}_{k})&\mathrm{if}\hskip 2.84544ptm\neq k\end{array}\right.$$

$$\displaystyle Cov(\bar{D}(\mbox{\boldmath$u$}_{m}),\bar{D}(\mbox{\boldmath$u$}_{k}))$$

$$\Vert {\hat{\Sigma }}_M-{\bar{\Sigma }}_M(K,{\sigma }^2){\Vert }_F^2=\Vert {\hat{\Sigma }}_M-\bar{Z}K{\bar{Z}}^{\prime }-{\sigma }^2\bar{V}{\Vert }_F^2$$

$$\hat{K}({\sigma }^2)=R^{-1}{Q}^{\prime }\left({\hat{\Sigma }}_M-{\sigma }^2\bar{V}\right)Q(R^{-1}{)}^{\prime }$$

$$\sum _{m,k}{\left(({\hat{\Sigma }}_M-Q{Q}^{\prime }{\hat{\Sigma }}_{M}Q{Q}^{\prime }{)}_{mk}-{\sigma }^2(\bar{V}-Q{Q}^{\prime }\bar{V}Q{Q}^{\prime }{)}_{mk}\right)}^2$$

$$\displaystyle F\succ 0\Leftrightarrow b<\frac{\mbox{\boldmath$e$}_{1}^{\prime}C\mbox{\boldmath$e$}_{1}}{\mbox{\boldmath$e$}_{1}^{\prime}D\mbox{\boldmath$e$}_{1}}$$

$$\displaystyle\sigma^{2}<\frac{\mbox{\boldmath$e$}_{1}^{\prime}R^{-1}Q^{\prime}\widehat{\Sigma}_{M}Q(R^{-1})^{\prime}\mbox{\boldmath$e$}_{1}}{\mbox{\boldmath$e$}_{1}^{\prime}R^{-1}Q^{\prime}\bar{V}Q(R^{-1})^{\prime}\mbox{\boldmath$e$}_{1}}$$

$$\displaystyle\hat{a}_{m}\approx 2^{-1/2}(\mbox{\boldmath$w$}^{\prime}_{m}\mbox{\boldmath$1$}_{n})^{1/2}/V_{D}(\mbox{\boldmath$u$}_{m})$$

$${\hat{K}}_a=R_a^{-1}{Q}_a^{{}^{\prime }}\left({\hat{\Sigma }}_M-{\sigma }^2\bar{V}\right)Q_a(R_a^{-1}{)}^{\prime }$$

$$\displaystyle\sigma^{2}_{a}<\frac{\mbox{\boldmath$e$}_{1}^{\prime}R^{-1}_{a}Q^{{}^{\prime}}_{a}\bar{A}^{1/2}\widehat{\Sigma}_{M}\bar{A}^{1/2}Q_{a}(R_{a}^{-1})^{\prime}\mbox{\boldmath$e$}_{1}}{\mbox{\boldmath$e$}_{1}^{\prime}R_{a}^{-1}Q_{a}^{{}^{\prime}}\bar{A}^{1/2}\bar{V}\bar{A}^{1/2}Q_{a}(R_{a}^{-1})^{\prime}\mbox{\boldmath$e$}_{1}}$$

$$\displaystyle\widehat{MSPE}(\Upsilon_{m})=\frac{1}{n_{t}}\displaystyle\sum_{i^{\prime}=1}^{n_{t}}(Y(\mbox{\boldmath$s$}_{i^{\prime}}^{*})-\hat{H}_{m}(\mbox{\boldmath$s$}_{i^{\prime}}^{*}))^{2}$$

$$F$$

$$E^{\prime }FE$$

$$\displaystyle\mbox{\boldmath$e$}_{1}^{\prime}F\mbox{\boldmath$e$}_{1}$$

$$\displaystyle\mbox{\boldmath$e$}_{1}^{{}^{\prime}}F\mbox{\boldmath$e$}_{1}$$

$$\displaystyle\mbox{\boldmath$e$}_{1}^{\prime}(C-bD)\mbox{\boldmath$e$}_{1}$$

$$b$$

$$\displaystyle b=\frac{v}{a(\mbox{\boldmath$e$}_{1}^{\prime}D\mbox{\boldmath$e$}_{1})}$$

$${\lambda }_1$$

$${\hat{\sigma }}_{g-1}^2={\hat{\sigma }}_g^2+{\delta }_g$$

$${\hat{K}}_g$$

$$\displaystyle\mbox{\boldmath$e$}^{*^{\prime}}_{g-1}\hat{K}_{g-1}\mbox{\boldmath$e$}^{*}_{g-1}$$

$$=$$

Full text

We present a result that shows when a matrix subtraction of a given form will give a positive definite matrix. Motivated by this result, we present an iterative Fixed Rank Kriging algorithm that ensures positive definiteness of the matrix required for prediction and show that under mild conditions the algorithm numerically converges. The modified Fixed Rank Kriging procedure is implemented to predict missing chlorophyll observations for very large regions of ocean color. Predictions are compared to those made by other well known methods of spatial prediction.

I Introduction

Gaussian process prediction (GPP) , require computations of $O(n^{3})$, which represents a formidable amount of computation when $n$ is very large. One way to overcome the computational burden has been using regression splines or penalized regression splines when the spatial dependence is modelled deterministically [1, 2]. In the case of the input dependence being modeled stochastically, several models have been suggested to reduce computational burden. Among them we have, tapering the covariance matrix [3], approximating the field using a Gaussian Markov random field [4], and doing computations in the spectral domain [5]. Often these stochastic models assume stationarity, an assumption that many times is justified for small sample sizes but likely to be too restrictive for large datasets. Moreover, published application of these methods primarily consider scenarios with a few thousand monitoring locations. When tens of thousands of observations or more are available, many of the techniques become unfeasible. In this article, we aim to modify a recent approach called Fixed Rank Kriging [6]. Particularly, we give a necessary and sufficient condition that ensures $F=C-bD$ is positive definite. Motivated by this result, we propose a new algorithm applied to Fixed Rank Kriging to estimate the necessary parameters $\sigma^{2}$ and matrix $K$ while ensuring positive definiteness of $\hat{K}$, the estimate of $K$. In this article we denote that a matrix $F$ is positive definite by $F\succ 0$. In sections II to II-B Fixed Rank Kriging is presented. This method makes predictions feasible even in the case of very large datasets, without assuming stationarity. In section III we introduce the new algorithm that ensures $\hat{K}\succ 0$ and we show that this algorithm converges numerically under general conditions. This is followed by an application using ocean color data.

II Fixed Rank Kriging

Cressie and Johannesson [6] introduced a new way of obtaining spatial predictions. This method is called Fixed Rank Kriging (FRK), and mainly addresses two drawbacks of classical GPP. One is that GPP cannot be implemented in the case of very large datasets since, in general, the computations involving the inverse of the covariance matrix, $\Sigma^{-1}_{K}$, are of the order $O(n^{3})$. Even a few thousand observations will make the GPP computations unfeasible. With FRK, the computations per prediction location are significantly faster while no longer requiring a stationary assumption. This section briefly reviews the inner workings of FRK as presented in [6], which from now on is referenced as CJ08. The model expressed in terms of a hidden process is,

[TABLE]

where H(\mbox{\boldmaths})=\mbox{\boldmathx}(\mbox{\boldmaths})^{\prime}\mbox{\boldmath\beta}+W(\mbox{\boldmaths}) is the hidden process of interest. \mbox{\boldmathx}(\mbox{\boldmaths})^{\prime}\mbox{\boldmath\beta} represents a random field mean and/or effects of $p$ covariates. For observations at $n$ locations \mbox{\boldmaths}_{1}),....,\mbox{\boldmaths}_{n}, let \mbox{\boldmathY}=(Y(\mbox{\boldmaths}_{1}),....,Y(\mbox{\boldmaths}_{n}))^{{}^{\prime}}, and $X$ be a $n\times p$ matrix with $k^{th}$ column (X_{k}(\mbox{\boldmaths}_{1}),....,X_{k}(\mbox{\boldmaths}_{n}))^{{}^{\prime}} for $k=1,..,p$. Assume $rank(X)=p$ with $p<n$. Furthermore, assume W(\mbox{\boldmaths}) and \epsilon(\mbox{\boldmaths}^{*}) are uncorrelated for locations \mbox{\boldmaths},\mbox{\boldmaths}^{*}\in\mathcal{D}\subset\Re^{d}. Also assume Var(\epsilon(\mbox{\boldmaths}_{i}))=\sigma^{2}v(\mbox{\boldmaths}_{i}) for $i=1,...,n$ and that W(\mbox{\boldmaths}) is a Gaussian process with Var(W(\mbox{\boldmaths}))<\infty \forall\mbox{\boldmaths}\in\mathcal{D}. If a fixed number $r$ of basis functions are chosen such that \mbox{\boldmathZ}(\mbox{\boldmaths})\equiv(Z_{1}(\mbox{\boldmaths}),...,Z_{r}(\mbox{\boldmaths}))^{\prime} are the basis functions evaluated at location $s$, CJ08 represent the spatial covariance of W(\mbox{\boldmaths}) between locations \mbox{\boldmaths}_{i},\mbox{\boldmaths}_{j} by,

[TABLE]

using a symmetric positive definite $r\times r$ matrix $K$, a matrix that will be estimated later on. Generally, (2) is not a function of distance between locations and is therefore a nonstationary covariance. CJ08 use the eigen-decomposition of $K$ to show that (2) can be interpreted as being similar to a truncated Karhunen-Loeve expansion but with non-orthogonal functions, \mbox{\boldmathZ}(\cdot).

Define $Z$ as the $n\times r$ matrix with $i^{th}$ row given by \mbox{\boldmathZ}(\mbox{\boldmaths}_{i})^{\prime}=(Z_{1}(\mbox{\boldmaths}_{i}),...,Z_{r}(\mbox{\boldmaths}_{i})), and let $V$ be a $n\times n$ diagonal matrix with v(\mbox{\boldmaths}_{i}) as the $i^{th}$ diagonal entry (v(\mbox{\boldmaths}_{i}) is assumed known). Combining (1) and (2) gives the covariance matrix of $Y$ in terms of $K$ as,

[TABLE]

By representing the $n\times n$ covariance matrix $\Sigma_{K}$ as in (3), the Sherman-Morrison-Woodbury equation may be used [7, page 50],

[TABLE]

Now calculation of $\Sigma_{K}^{-1}$ requires only the inversion of a $r\times r$ matrix $K$ and a $n\times n$ diagonal matrix $V$. Predictions at locations \mbox{\boldmaths}_{o} are now possible by using GPP with $\Sigma_{K}^{-1}$ specified as in (4). Thus by substitution into the GPP equation accounting for measurement error [8] we get,

[TABLE]

where \hat{\mbox{\boldmath\beta}}=(X^{{}^{\prime}}\Sigma_{K}^{-1}X)^{-1}X^{{}^{\prime}}\Sigma_{K}^{-1}\mbox{\boldmathY}, while the mean squared prediction error of \hat{H}(\mbox{\boldmaths}_{o}) is,

[TABLE]

Both the covariance parameter matrix $K$ and the scalar $\sigma^{2}$ need to be estimated to proceed with FRK.

II-A Basis functions

Smoothing spline bases, regression spline bases, radial basis functions, eigenvectors or wavelets are among some basis functions that could be used with FRK. In this work, we use the same bisquare basis functions as [6]. These revolve around the concept of gridding the locations at different resolutions and obtaining centroid locations. Basis functions at a coarse resolution/scale, capture general global attributes of the process. The finer the scale becomes, the more local are the attributes captured by each basis function at that resolution. Another benefit of the bisquare basis function is that is has an intuitive interpretation, the closer a location is to a centroid, then the closer to 1 is the basis function, while the further they are, the closer to zero is the basis function. Furthermore each basis function has local support. For prediction, computations of $Z^{\prime}A^{-1}Z$ and Z^{\prime}\mbox{\boldmatha} for any $A^{-1}$ and $a$ are $O(nr^{2})$ for any basis function where $r$ is the total number of basis functions [6]. But when matrix $Z$ is sparse, as is the $Z$ matrix resulting from using bisquare basis functions, then the operations required for prediction are $O(kr^{2})$ for $k<r<<n$.

Suppose a FRK model is designed for a square region of $3,600$ regularly spaced locations such that the covariance (2) is composed of two scales of variation, namely a coarse scale with 4 functions, and a finer scale with 25 functions. Figure 1 displays the basis functions at the coarse scale. Similarly, Figure 2 displays the basis functions at the fine scale.

II-B Parameter estimation

Positive definiteness of the covariance matrix is required to ensure invertibility of $\Sigma_{K}^{-1}$ and positive mean squared prediction errors [8]. CJ08 bin the data for the purpose of estimation of the covariance. This way, the estimation and fitting of the covariance will depend on the number of bins $M$, and not on the number of locations $n$, where $M<<n$. Assigning bin centers \{\mbox{\boldmathu}_{m}:m=1,...,M\}, a neighborhood N(\mbox{\boldmathu}_{m}) is defined. The neighborhood N(\mbox{\boldmathu}_{m}) could be set up by some Euclidean distance as a threshold (independent of direction, a sort of circular neighborhood, or by distances in the horizontal and vertical direction). Then define the neighborhood indicator variable as,

[TABLE]

Note, that a missing value at location \mbox{\boldmaths}_{i} would imply a weight of zero. Furthermore, designate \mbox{\boldmathw}_{m}=(w_{m1},...,w_{mn})^{{}^{\prime}}, and let \mbox{\boldmath1}_{n} be the $n$ dimensional column vector of ones. Without any knowledge of the process we may use the ordinary least squares estimator \hat{\mbox{\boldmath\beta}} of $\beta$, to detrend the data such that, D(\mbox{\boldmaths}_{i})=Y(\mbox{\boldmaths}_{i})-\mbox{\boldmathx}(\mbox{\boldmaths}_{i})^{{}^{\prime}}\mbox{\boldmath\hat{\beta}}. Let \mbox{\boldmathD}=(D(\mbox{\boldmaths}_{1}),...,D(\mbox{\boldmaths}_{n}))^{{}^{\prime}}.

Then the average of the residuals within the $m^{th}$ bin will be,

[TABLE]

Let \mbox{\boldmath\bar{D}}=(\bar{D}(\mbox{\boldmathu}_{1}),...,\bar{D}(\mbox{\boldmathu}_{M}))^{{}^{\prime}}. An estimate of the true covariance matrix of the binned residuals, \Sigma_{M}=Var(\mbox{\boldmath\bar{D}}) based on the empirical method of moments estimator has element $(m,k)$ given by,

[TABLE]

for $m,k\in\{1,...,M\}$ where V_{D}(\mbox{\boldmathu}_{m})=\sum_{i=1}^{n}w_{mi}D(\mbox{\boldmaths}_{i})^{2}/\mbox{\boldmathw}_{m}^{{}^{\prime}}\mbox{\boldmath1}_{n} and C_{D}(\mbox{\boldmathu}_{m},\mbox{\boldmathu}_{k})=\bar{D}(\mbox{\boldmathu}_{m})\bar{D}(\mbox{\boldmathu}_{k}). Matrix $\widehat{\Sigma}_{M}$ has the $(m,k)$ elements seen in equation (12).

The intention is to approximate $\hat{\Sigma}_{M}$ by a function of spatial dependence parameters $K$ and $\sigma^{2}$. That is, we want a matrix $\bar{\Sigma}(K,\sigma^{2})$ with the following $\bar{\Sigma}(K,\sigma^{2})_{mk}$ entries,

[TABLE]

such that if $\bar{Z}=(\bar{Z}^{{}^{\prime}}_{1},....,\bar{Z}^{{}^{\prime}}_{M})^{{}^{\prime}}$ is the $M\times r$ matrix of binned basis functions, and $\bar{V}=diag(\bar{V}_{11}^{*},...,\bar{V}_{MM}^{*})$, then $\bar{\Sigma}(K,\sigma^{2})=\bar{Z}K\bar{Z}^{{}^{\prime}}+\sigma^{2}\bar{V}$222Instead of deriving $\bar{V}_{mk}^{*}$ strictly from (II-B), we use \bar{V}_{mk}=\mbox{\boldmathw}_{m}^{{}^{\prime}}V\mbox{\boldmathw}_{k}/(\mbox{\boldmathw}_{m}^{{}^{\prime}}\mbox{\boldmath1}_{n}) to get a more appropriate estimator of $\sigma^{2}$. Although the $\sigma^{2}$ estimator has been rescaled, the resulting $\hat{K}$ remains invariant to the rescaling. As a consequence, FRK predictions and their standard errors are unaffected.. In order to approximate $\hat{\Sigma}_{M}$ by $\bar{\Sigma}(K,\sigma^{2})$, CJ08 quantify the ’best’ approximation in terms of the Frobenius norm. Specifically for FRK,

[TABLE]

is minimized to obtain an estimate of $K$ and $\sigma^{2}$. To ease matrix computations, the QR-decomposition of $\bar{Z}$ is performed. That is, $\bar{Z}=QR$ where $Q$ is a column-wise orthonormal matrix of dimension $M\times r$, and $R$ is a nonsingular $r\times r$ upper triangular matrix. Then the $K$ that minimizes (14) is given by,

[TABLE]

Therefore $\hat{K}(\sigma^{2})$ depends on $\sigma^{2}$. Substituting (15) into $\bar{\Sigma}_{M}(K,\sigma^{2})$ alters (14) such that $\hat{\sigma}^{2}$, is the estimator of $\hat{\sigma}^{2}$ that minimizes,

[TABLE]

Note that remarkably, the minimization problem for $\sigma^{2}$ is clearly of the form of a simple least squares problem with slope $\sigma^{2}$ and no intercept and is therefore easy to obtain (of course $\hat{\sigma}^{2}$ must be constrained to be positive since $\sigma^{2}$ is a variance parameter)333Although [6, page 216] state that $r\leq M<n$, for the estimation of the parameters we must have $r\neq M$. If the number of bins $M$ were equal to the number of basis functions $r$ this would imply that also $QQ^{\prime}=I$ which would make the minimization problem (14) impossible to solve. Therefore for FRK, it is required that $r<M<n$ ($r<M$ to ensure column-wise orthogonality of $Q$).

When $\sigma^{2}=0$, $\hat{K}$ is positive definite while if $\sigma^{2}>0$ the matrix is positive definite conditional on the value of $\sigma^{2}$. There is no guarantee that the estimate $\hat{\sigma}^{2}$ will result in a positive definite matrix. In fact, Shi and Cressie (2007, page 671) state that sometimes $\hat{\sigma}^{2}$ needs to be adjusted so that $\hat{K}$ is positive definite but do not give further details. In this work, an algorithm is proposed to use the positive definiteness of $\hat{K}$ as a constraint.

We simulate a Gaussian mean [math] stationary Gaussian process with an exponential covariance function, a partial sill of 5.5, range of 1. The nugget was chosen to be 1.375, so that the signal to noise ration was 4. A sequence of 1,500 potential values of $\sigma^{2}$ were used to then calculate $\hat{K}(\sigma^{2})$, its smallest eigenvalue $\lambda_{min}$, and the sum of squares produced in (16) by the $\sigma^{2}$ value. $\sigma^{2}=2.4307$ is the result of minimizing (16) when only constrained so that the estimator is positive; whereas $\sigma^{2}=1.6856$ gives the smallest sum of squares such that $\lambda_{min}>0$. Figure 3 shows how $\lambda_{min}$ and the sum of squares (16) related to $\sigma^{2}$. The issue at hand is clearly seen, the minimization of the sum of squares to find $\hat{\sigma}^{2}$ requires the constraint that $\hat{K}$ is positive definite. Hence the solution is a subspace of the unconstrained sum of squares problem for $\hat{\sigma}^{2}$. Yet, positive definiteness is a quadratic condition and appears tricky to implement as a constraint. However, it turns out that the positive definiteness of $\hat{K}$ can be applied as a upper bound constraint on $\hat{\sigma}^{2}$.

III A new algorithm to estimate the FRK parameters

Clearly from (15), if $\hat{\sigma}^{2}$ is too large, then $\hat{K}$ will have negative diagonal entries. A symmetric matrix with negative entries in the diagonal cannot be positive definite [7, page 141]. We wish to impose positive definiteness as a constraint when minimizing (16). The following lemma leads us in the right direction,

Lemma 1**.**

*Define $F=C-bD$ where all matrices $F,C,D$ are $r\times r$ real matrices, $C\succ 0$, and $D\succ 0$, and $C$ and $D$ are symmetric. Furthermore, assume $F$ has distinct eigenvalues and that $b$ is any constant such that $b>0$. Then,

[TABLE]

where \mbox{\boldmathe}_{1} is the eigenvector associated with the minimum eigenvalue of $F$, $\lambda_{1}$.

The interpretation of the lemma is that, a subtraction of two positive definite symmetric matrices of a certain form gives us a matrix $F$ that is positive definite under certain conditions, among them that the scalar $b$ is not too big. The following Corollary shows how $\hat{K}$ is a special case of Lemma 1.

Corollary 1**.**

*Assume $\hat{K}$ in (15) has distinct eigenvalues, $\lambda_{1}<....<\lambda_{r}$. Then $\hat{K}$ is positive definite if and only if,

[TABLE]

where \mbox{\boldmathe}_{1} is the $r\times 1$ normalized eigenvector corresponding to the smallest eigenvalue $\lambda_{1}$ of $\hat{K}$.

Result (17) motivates the use of the positive definiteness requirement of $\hat{K}$ as a linear constraint on $\hat{\sigma}^{2}$. The following algorithm iteratively estimates $\sigma^{2}$ and $K$.

**FRK parameter estimation algorithm

**

Calculate $Q$, $R$, $\bar{V}$, and $\hat{\Sigma}_{M}$ as described in section II-B. 2. 2.

Estimate $\sigma^{2}$ by minimizing equation (16) only subject to a constraint that $\hat{\sigma}^{2}>0$. Start at zero an index of the iteration, $g=0,1,..$. Set $\hat{\sigma}^{2}_{g}$ as the result from (16). 3. 3.

Calculate $\hat{K}_{g}\equiv\hat{K}(\hat{\sigma}_{g})$ using (15). 4. 4.

Check if $\hat{K}_{g}\succ 0$. This is so if $\lambda_{min,g}>0$. If it is, we stop here. If it is not, calculate an upper bound according to (17) for $\hat{\sigma}_{g}^{2}$. Let the upper bound be $\hat{\sigma}_{u,g}^{2}$. 5. 5.

Minimize (16) over $\hat{\sigma}^{2}_{g}$ but now subject to both, the greater than zero constraint and to the upper bound $\hat{\sigma}_{u,g}^{2}$ constraint. 6. 6.

Repeat steps 3-5 above until $\hat{K}_{g}\succ 0$.

This algorithm is a special case of the cutting plane algorithm developed by Tuy [9, 10]. What remains is to show that as $g^{th}$ increases, $\lambda_{min,g}$ of $\hat{K}_{g}$ increases while the upper bound provided by (17) decreases.

Now we are ready to show the numerical convergence of the proposed FRK algorithm by stating the following theorem,

Theorem 1**.**

If $\lambda_{min,g}$ is the minimum eigenvalue of $\hat{K}_{g}$ at iteration $g$, $\hat{K}_{g}$ has distinct eigenvalues $\lambda_{min,g},...,\lambda_{max,g},\forall g$ and $\sigma_{u,g}^{2}$ is the upper bound found in Step 4 of the FRK parameter estimation algorithm at iteration $g$, then $\lambda_{min,g}>\lambda_{min,g-1}$ if and only if $\hat{\sigma}_{g}^{2}<\hat{\sigma}_{g-1}^{2}$.

Returning to the mean [math] stationary Gaussian process with an exponential covariance function adopted earlier, implementing our FRK parameter estimation algorithm using $M=100$ results in $\hat{\sigma}^{2}=1.6861$ after 8 iterations. Figure 4 demonstrates how $\lambda_{min,g}$ and the sum of squares from (16) increase per iteration $g$, while the number of negative eigenvalues and $\hat{\sigma}_{g}^{2}$ decrease with every iteration.

Kang and Cressie [11] propose an empirical way to ensure that $\hat{K}$ is positive definite by updating the eigenvalues of $\widehat{\Sigma}_{M}$. To maintain the variability features of the original estimate $\widehat{\Sigma}_{M}$ they must choose specific values of a parameter in their algorithm. In our algorithm we do not change the estimate of spatial covariance and instead iterate among the estimates of $\hat{\sigma}^{2}$ and $\hat{K}$. Now that we can ensure that $\hat{K}_{g}\succ 0$ the following steps are taken to attain the FRK prediction of a random field of interest:

Calculate the ordinary least square estimate \mbox{\boldmath\hat{\beta}}=(X^{{}^{\prime}}X)^{-1}X^{{}^{\prime}}\mbox{\boldmathY} of $\beta$. 2. 2.

Obtain residuals \mbox{\boldmathD}=\mbox{\boldmathY}-X\mbox{\boldmath\hat{\beta}} according to the ordinary least squares estimate of $\beta$. 3. 3.

Set up the matrix of basis functions $Z$. 4. 4.

Next the stage is set to estimate the covariance matrix $\Sigma_{K}$ as in (3). 5. 5.

Calculate the method of moments estimate $\widehat{\Sigma}_{M}$ of $\Sigma_{M}$ using \bar{D}(\mbox{\boldmathu}_{m}) 6. 6.

Obtain $\bar{Z}$ and $\bar{V}$. Then perform the QR-decomposition of $\bar{Z}$. 7. 7.

Estimate $\sigma^{2}$ and $\hat{K}$ according to the algorithm proposed in this section and shown to numerically converge in Theorem 1. 8. 8.

Predict the random field H(\mbox{\boldmaths}_{o}) using (5) and its standard error with the square root of (6).

Several things are worth mentioning. First, it is implicitly assumed that the binning has been performed in such a way that missing values are not present in the binned version of the data. Otherwise weighting of the covariance estimates would not be possible (an empty bin results in V_{D}(\mbox{\boldmathu}_{m})=\infty for that bin as seen from (12)). More importantly in the presence of missing values the positive definiteness of the matrix could no longer be guaranteed. A remedy could be to remove empty bins from (12) or impute them by some rule. $x,y$ coordinate based averages or local averaging are just some possibilities of imputing missing values. Cressie [8] discusses median polish, another tool that could be used for imputation of missing bins. Little and Rubin [12] provide a much broader emphasis on ways to deal with missing values.

Another noteworthy fact is that storage of the $n\times n$ $\Sigma^{-1}_{K}$ is often required in classical geostatistics, a task of the order of $O(n^{2})$. In addition to savings on computational time, FRK provides substantial reduction in storage.

Another point worth mentioning is that sometimes some bins may be considerably more variable than others. To account for differences in variability among the bins, the influence of potential outliers, and for bins that have more data, the parameter estimation can be weighted. This is analogous to the use of weighted least squares instead of ordinary least squares in variogram fitting for the same purpose. CJ08 also took this into consideration and based on [see 8, page 95] and references therein suggest,

[TABLE]

as a weight. This translates into a rescaling of the minimization problem (14). If $\bar{A}=diag(\hat{a}_{1},....,\hat{a}_{M})$, then the norm is equivalent to $\|\bar{A}^{1/2}\hat{\Sigma}_{M}\bar{A}^{1/2}-\bar{A}^{1/2}\bar{\Sigma}_{M}(K,\sigma^{2})\bar{A}^{1/2}\|_{F}^{2}$ and,

[TABLE]

where we used $\bar{A}^{1/2}\bar{Z}=Q_{a}R_{a}$. Therefore virtually no computational cost is added by using the weights. Furthermore, the previous results in this article still hold.

Corollary 2**.**

*Assume $\hat{K}_{a}$ in (18) has distinct eigenvalues, $\lambda_{1}<....<\lambda_{r}$. Then $\hat{K}_{a}$ is positive definite if and only if,

[TABLE]

where \mbox{\boldmathe}_{1} is the $r\times 1$ normalized eigenvector corresponding to the smallest eigenvalue $\lambda_{1}$ of $\hat{K}_{a}$.

The proof is almost identical to Lemma 1 and Corollary 1. Theorem 1 also still holds.

Fixed Rank Kriging is a method intended for very large datasets. Comparisons with GPP are made difficult for two reasons. First, if the dataset is too large, classical GPP is not feasible. On the other hand if a dataset is small enough for GPP, it might be too small to estimate the spatial covariance properly using (12). For example if $n=3,600$ one should probably have at least 30 pixels per bin, which doesn’t leave much room for a multiresolution basis.

Finally since the inverse computation is not an $O(n^{3})$ anymore when FRK is used, a Likelihood approach that would more fully take into account the probability distribution of the data in comparison to the method of moments approach is plausible. Yet keep in mind that in this case one of the parameters of interest, $K$, is a matrix. Stein [13] fits a covariance of the form (3) by parameterizing $K$ and then maximizing the likelihood.

IV Chlorophyll data

Ocean color observations enable scientists to study several biological and biogeochemical properties of the oceans. In part, ocean color can measure surface phytoplankton (microscopic ocean plants) since color in most world’s oceans in the visible light region (wavelength of 400-700nm) varies with the concentration of chlorophyll and other components, i.e the more phytoplankton present, the greater the concentration of plant pigments, ergo the greener the water [14]. Ocean color is crucial for: the study of organic matter produced by algae and bacteria, the study of the biochemistry of the ocean, the assessment of the role of the ocean in the carbon cycle, and the potential global warming trend.

Based on prior and ongoing ocean color satellite missions, scientists can now study the spatial and temporal variability of the biological, chemical and physical processes that regulate ocean color around the globe. For example, [15] and [16] present studies of the spatial correlation of chlorophyll at the mesoscale (about 10-200km). Substantial progress has been made in analyzing these ocean color satellite data, See [17, 14, 18]. But many computationally intensive statistical challenges remain. The ocean color satellite datasets are massive, with hundreds of thousands of observations or more around the globe. Ocean processes are generally non-stationary in both space and time. Furthermore, due to many factors, satellite data comes with large amounts of missing data.

In this section our goal is to compare the prediction capability of of the modified Fixed Rank Kriging (FRK) algorithm , to other well known prediction methods, namely: ordinary least squares, regression thin plate splines, and universal kriging, whenever computations required for the latter are feasible.

IV-A Satellite data

Satellites allow us to gather environmental data on a truly global scale. Currently, there are several satellites sending ocean color data to different agencies from several countries (See http://www.ioccg.org/sensors_ioccg.html). NASA is in charge of two of these satellites that carry instruments sending ocean color data. The Sea-viewing Wide Field-of-view Sensor (SeaWiFS) is one of them, and was deployed in 1997. The other is the Moderate Resolution Imaging Spectroradiometer (MODIS-Aqua), an instrument available in the Aqua satellite operating since 2002. These sensors offer different designs, orbits, and accuracy of ocean color data. Specifically, the sensors measure the water leaving radiance, $L_{wN}(\lambda)$, a subsurface radiance reflected out of the ocean through an air-sea interface. Both SeaWiFS and MODIS-Aqua measure the water leaving radiance $L_{wN}(\lambda)$ at several wavelengths $\lambda$. These two instruments offer different space-time resolutions, different sampling patterns, and different measurement uncertainties. Table I summarizes some traits that characterize each of the two sensors.

The satellite raw data plus instrument data (’Level 1’) is calibrated and reprocessed to give geophysical values, Level 2, derived after applying several algorithmic and atmospheric adjustments (http://oceancolor.gsfc.nasa.gov). The Level 2 data are spatially and temporally averaged to give Level 3 statistical data. For example, SeaWiFS has a $9km\times 9km$ spatial resolution at the equator for Level 3 products while MODIS has $4.6km\times 4.6km$ spatial resolution at the equator for Level 3 products. See [20].

The $L_{wN}(\lambda)$ data used for our analysis is processed by a semi-analytical model developed by Garver, Siegel and Maritorena, GSM01 [18]. The algorithm expresses the water leaving radiance in terms of CHL and the Inherent Optical Properties (IOP): CDM and BBP. In short, the GSM01 model inverts observations of the normalized water leaving radiance spectrum, $L_{wN}(\lambda)$, to estimate, Chl, CDM and BBP [21].

Data from two NASA satellites are available. Specifically for SeaWiFS, the time period of September 04, 1997 through July 04, 2007 is available. Aqua results from July 04, 2002 through July 04, 2007 are also available. Missing observations result from orbital sampling, sun glint and cloud cover. Regions of the globe often have observations for only about 30% of the days (sometimes less) across all years. Figure 5 shows the satellite measurement pattern for the Aqua sensor for one day (December 31, 2002), and highlights the space-time missing data problem. The orbital track of the satellite is clearly visible (especially in the Southern Hemisphere), yet on occasion, regions that fall along the sampling track do not have observations. In fact, there are very few observations north of $40^{o}$N, while more data is present south of $40^{o}$S. This pattern in the missing values is due to a seasonal effect on the mechanism of missing values. We focus on open ocean waters where GSM01 output is considered to be more reliable than that coming from coastal water reflectance [21].

Maritorena and Siegel [22] combine the $L_{wN}(\lambda)$ values from SeaWiFS and Aqua combined to produce CHL and IOP data. This results in better spatial/temporal coverage and sometimes lower uncertainty surrounding the variables obtained (pixels with multiple observations are imputed a weighted average of satellite data according to their respective uncertainties). More recently [23] included MERIS data as well in the merging procedure. We begin the analysis on the field of chlorophyll values and the chlorophyll dependence on spatial location.

V Chlorophyll for a very large region in the North and

South Pacific

In this section, we examine the predictions obtained from four prediction methods: ordinary least squares (OLS), an additive model (AM), and fixed rank kriging (FRK). For this analysis, the observations consist of 8-day Sea-viewing Wide Field-of-view Sensor (SeaWiFS) and Moderate Resolution Imaging Spectroradiometer (Modis) Aqua merged data resulting from the GSM01 algorithm. Campbell [24] discusses how chlorophyll follows approximately a lognormal distribution. We follow this result and model the chlorophyll on the natural logarithm scale and denote Y(\mbox{\boldmaths}_{i})=log(CHL_{i}) for location (pixel) $i$. Hypothesis tests on $\beta$ in a preliminary OLS model with \mbox{\boldmathx(s_{i})}=(1,s_{i1},s_{i2})^{\prime} as the $i^{th}$ row of the $n\times 3$ matrix $X$ where \mbox{\boldmaths}_{i}=(s_{i1},s_{i2}) for $i=1,\ldots,n$ (where $n=3600$), and with \mbox{\boldmath\beta}=(\beta_{1},\beta_{2},\beta_{3})^{\prime}, implied a significant effect in both North-South and East-West directions.

Due to the size of the data, the AM is implemented using thin plate regression splines. More precisely, the basis functions for the thin plate regression splines model are chosen by eigen-truncating the matrix with entries coming from the basis functions of thin plate splines [2]. The FRK procedure is implemented in Matlab using a 32-bit Linux machine.

V-A Measure of predictor performance

$n_{t}$ observations are put aside for testing the performance of a predictor. The remaining observations are used to fit OLS, AM, and FRK models. Then the estimated mean squared prediction error can be estimated by cross validation,

[TABLE]

calculated at test locations \mbox{\boldmaths}_{i^{\prime}}^{*} for $n_{t}$ prediction locations, and $\Upsilon_{m}$ is the $m^{th}$ modeling procedure. We also note that (19) must be used with caution. Care is needed, since this method does not prove that a spatial model is correct, only that it is ’not grossly incorrect’ [8].

V-B Satellite data for regions of interest

The regions in the analysis are 130-155W by 5-30N in the North Pacific and 125-150W by 5-30S in the South Pacific, both which include 90,000 locations. Kriging would require the storage of a $90000\times 90000$ covariance matrix for each region and the inversion of this matrix is $O(90000^{3})$. In the North Pacific region the 8-day data has 1947 missing values for the time period starting at Julian day 73 of year 2003, and in the South Pacific 2,403 missing values for the time period starting at Julian day 177 of year 2007. These time periods are conveniently chosen so that we can separate different sizes of test data for the comparison of the prediction models like we did in the previous section. Figure 6 displays the South Pacific region of interest. Specks of pixels with apparent high values relative to neighbors can be seen around $15^{o}$S, $145^{o}$W. These are due to the coasts of the French Polynesia islands in that region. Clearly most of the variability occurs in a North-South direction, with an increasing trend in the northward direction. The North Pacific region also shows most of the variability in the North-South direction (but with an increasing trend in the southward direction). Based on these preliminary images, we assume a linear trend with the coordinates in each direction.

The AM is implemented using thin plate regression splines. Exploratory analysis lead us to choose 100 basis functions to fit the thin plate regression spline model.

Fixed Rank Kriging is implemented using bisquare basis function with three scales of variation 16, 64, and 225 functions respectively. Parameters $K$ and $\sigma^{2}$ were estimated using $M=900$, $V=I$ and following the new procedure detailed in section II.

Missing values are omitted before the comparison of the predictors. Then 15% of the data are randomly removed to test each model. The remaining observations are used to fit OLS, AM and FRK models and then each model is used to predict the test data. The procedure is repeated 50 times. Then the $\widehat{MSPE}$ is obtained using (19). The study was conducted also using 25% and 50% of the observations as test data. Table II presents the results of the analysis for the North Pacific. While Table III does so for the South Pacific. OLS gives the highest estimates of $MSPE$. For both the North and South Pacific we find that FRK outperforms not only OLS but AM as well. The FRK method takes about 167 seconds to calculate $\hat{K}$ and $\hat{\sigma}^{2}$ using the procedure introduced in section III and estimate the process at all 44,026 locations in the North Pacific (see Table II) when 50% of the original data is used to test the prediction method. A regression thin plate spline implements a basis function smoothed equally in all directions. Anisotropy in the large region could hinder any method based on such an assumption. Yet a tensor product of regression splines [2], in both spatial directions did not improve on the results of the thin plate regression splines shown here.

Moreover, the mean $\widehat{MSPE}$ of the prediction methods almost do not change as more missing values are present in the data. In the case of the OLS and AM this may be so because, as suspected from Figure 6, most of the variability occurs in the North-South direction and, if thought of deterministically, the association between chlorophyll and location is almost linear. Even when half of the data is removed for cross-validation, the amount of remaining observations may be enough that the FRK predictor performs well. Further study would be needed to confirm this.

FRK allows quick prediction of missing values in a massive region without assuming that the spatial association is stationary. Figure 7 displays how the ocean color image would look like if missing values where predicted using FRK and imputed into the original image (right panel). The prediction method does not induce any unwanted discontinuities or extreme outliers in the image. Unfortunately, with the number of observations used in this section GPP cannot be used. We suspect that it will underperform FRK because in this scenario the spatial association may vary considerably across the region.

VI Final Remarks

FRK can be viewed as a Gaussian process random effects model that allows prediction in the case of very large datasets without assuming stationarity and without storing the $n\times n$ covariance matrix. We give a result showing when the difference of two matrices (of a certain form) is guaranteed to be a positive definite matrix. Based on that result, we have proposed an iterative algorithm that finds $\hat{\sigma}^{2}$ such that $\hat{K}$ is ensured to be positive definite. We prove that the algorithm always converges numerically. The results hold when parameter estimation is weighted according to bin variability.

In this work we implement the modified FRK algorithm to ocean color data. We find that the prediction method gives smooth predictions without assuming stationary spatial dependence and can make these predictions quickly, even when the data are very large. In comparison, thin plate regression splines require more computational time while, depending on the pattern and amount of missing data, not performing as well as FRK. Moreover, the FRK predictions appear to hold the association between CHL and the IOPs (results not shown).

Many directions are possible for future work. It would be useful to see the performance of FRK when different basis functions and different values of $r$ are used. Shi and Cressie [25] use wavelets as basis functions. Yet efficient implementation of wavelets requires the absence of missing values. The authors impute the missing values to obtain the basis function matrix. How this step affects their final results is unclear and further study is warranted. We conjecture that a flexible basis function can be derived with the use of truncated tensor product basis functions. Spatio-temporal and multivariate extensions would also be useful. Other explanatory variables known to affect ocean color could be included in the prediction model: sea surface temperature, wind and in situ measurements of ocean color could be used. The effect of some of the properties ($M$, $r$, basis function class, etc.) of the FRK model on ocean color predictions needs further study. The FRK predictions of missing values will provide more coverage of the ocean. The biological, chemical and physical forcings could be studied at different spatial scales, taking into account the predicted observations and the uncertainty surrounding their estimation.

Appendix A Proofs

A-A Proof of Lemma 1

Proof.

Since $F$ is symmetric (it is the result of subtracting two symmetric matrices), by the spectral theorem,

[TABLE]

where $\Lambda=diag(\lambda_{r},...,\lambda_{1})$ is the matrix with eigenvalues of $F$ and $E$ is the associated eigenvector matrix. Equality in (20) holds since the eigenvalues are all distinct, implying the eigenvector matrix is orthogonal, $E^{\prime}E=I\Rightarrow E^{\prime}=E^{-1}$. Note that (20) implies that:

[TABLE]

Assume $F\succ 0$ to prove the upper bound on $b$ and then complete the proof by assuming the upper bound in $b$ and showing that then $F\succ 0$. Now, given that $F$ is positive definite, then all its eigenvalues are positive. Therefore by (21),

[TABLE]

To complete the proof, assume b<\frac{\mbox{\boldmathe}_{1}^{\prime}C\mbox{\boldmathe}_{1}}{\mbox{\boldmathe}_{1}^{\prime}D\mbox{\boldmathe}_{1}}. Then $\exists\hskip 2.84544pta\in\Re:a>1$ such that,

[TABLE]

where v=\mbox{\boldmathe}_{1}^{\prime}C\mbox{\boldmathe}_{1}. Note that $v>0$ since $C\succ 0$. Then from (21),

[TABLE]

Substituting (22) for $b$ gives,

[TABLE]

The last inequality is due to $v>0$ and $a>1$. ∎

A-B Proof of Corollary 1

Since, as stated in section II, $\hat{\Sigma}_{M}$ and $\bar{V}$ are positive definite, $R^{-1}Q^{\prime}\widehat{\Sigma}_{M}Q(R^{-1})^{\prime}$ and $R^{-1}Q^{\prime}\bar{V}Q(R^{-1})^{\prime}$ are also positive definite [7, p. 141]. The proof of the corollary is exactly the same as for Lemma 1 with $F=\hat{K},C=R^{-1}Q^{\prime}\widehat{\Sigma}_{M}Q(R^{-1})^{\prime},D=R^{-1}Q^{\prime}\bar{V}Q(R^{-1})^{\prime}$, and $b=\hat{\sigma}^{2}$.

A-C Proof of Theorem 1

Proof.

Assume $\lambda_{min,g}>\lambda_{min,g-1}$ and denote the minimum normalized eigenvector as \mbox{\boldmathe}^{*}=\mbox{\boldmathe}_{min}, then the eigenvector associated with $\lambda_{min,g}$ at iteration $g$ will be \mbox{\boldmathe}^{*}_{g}. Furthermore, let $\hat{K}_{g}=C-\hat{\sigma}^{2}_{g}D$ where we define $C=R^{-1}Q^{\prime}\widehat{\Sigma}_{M}Q(R^{-1})^{\prime}$, and $D=R^{-1}Q^{\prime}\hat{V}Q(R^{-1})^{\prime}$. Showing $\hat{\sigma}^{2}_{g}<\hat{\sigma}^{2}_{g-1}$ is equivalent to demonstrating that if,

[TABLE]

then $\exists\hskip 2.84544pt\delta_{g}\hskip 2.84544pt\in\Re:\delta_{g}>0$ for all $g$. Using (23) one can write,

[TABLE]

Now one has everything required for the first part of the proof. Assuming that all eigenvalues of $\hat{K}_{g}$ are distinct, then,

[TABLE]

The first equality comes from (21) with $F=\hat{K}$. The inequality afterwards is by assumption in this part of the proof. The inequality in (25) is due to the Rayleigh-Ritz theorem result [26]. (A-C) is obtained by substituting $\hat{K}_{g}$ with (24). From the last equation we see that, given that $D$ is known to be positive definite, $\delta_{g}>0$ for all $g$ therefore showing that $\hat{\sigma}^{2}_{g}<\hat{\sigma}^{2}_{g-1}$ by (23).

To complete the proof, now assume $\hat{\sigma}^{2}_{g}<\hat{\sigma}^{2}_{g-1}$ and write $\hat{\sigma}^{2}_{g}=\hat{\sigma}^{2}_{g-1}-\delta_{g}$. From (15) we may write,

[TABLE]

Then,

[TABLE]

from this result we deduce that $\lambda_{min}(\hat{K}_{g})=\lambda_{min}(\hat{K}_{g-1}+\delta_{g}D)$. As a final step,

[TABLE]

the last inequality is due to theorem 8.1.5 in [7, p. 396]. ∎

A-D Proof of Corollary LABEL:cor:othercor and Corollary LABEL:theorem:iterationalgweighted

Proofs are almost identical to Lemma 1 and Theorem 1 respectively.