Imputing missing groundwater observations

The paper considers a network of observation wells which are used for monitoring the hydrologic situation in a certain aquifer. The water levels which are measured in the wells form the basis of regular reports, numerical modeling, water resources assessments, and various research projects. Many of these applications are based on annual records of data. That is, data records which are collected once a year, at about the same time. In Israel, for example, each year is composed of two seasons, a rainy one and a dry one, and the annual records of water levels are measured in October, at the end of the dry season. The term water level refers to the height of the water table at a well. Assume that annual records of water levels were taken during m years, in a group of n neighboring wells. Then the collected data form an m x n matrix, in which each row refers to a different year, and each column refers to a different well. Let A denote the resulting data matrix, and let a_ij denote the (i,j) entry of A. Then a_ij gives the water level measured at the ith year in the jth well, i = 1, …, m, j = 1, …, n (see Table 1). Most reports and models require a complete data matrix as input. That is, a matrix without missing entries. Thus, even one missing entry may prevent us from using the collected data. However, as explained below, in practice it is quite common that some entries of the matrix are missing. In this case it is desired to assign ‘appropriate values’ to the missing entries. Once all the missing water levels are imputed, the data matrix can be used for any purpose.

The task of imputing missing values is, therefore, to design an algorithm that substitutes appropriate values into the missing entries of a given m x n data matrix, A = (a_ij). Other common names of this task are ‘missing values estimation’, ‘completing missing entries’, ‘matrix completion’, ‘recovering’, ‘reconstructing’, etc. The need for completing missing entries of a data matrix is a well known problem that arises in many applications. Many ingenious algorithms have been proposed and there is vast literature on this issue. In particular, the statistical treatment of missing data is a relatively matured area which preceded most of the other fields, e.g., Afifi & Elashoff (1966) or Dempster et al. (1977). Other so-called ‘traditional’ fields include business reports, psychometrica, biology and meteorology. Modern applications arise in machine learning, data mining, DNA microarray data, computer vision, recommender systems, and collaborative filtering. See Dax (2014) for a recent comprehensive survey of imputing methods. Most of the methods consider the imputing problem within the context of a specific application. The current research concentrates on imputing methods which are able to take advantage of the special features that characterize annual records of water levels.

There are several reasons for having missing records of water levels. One comes from human factors such as shortage of manpower, illness, strikes, typing errors, etc. A second stems from (unexpected) technical difficulties, such as technical faults in the measuring equipment, blocked roads, cars breakdown, etc. Another difficulty comes from the fact that water levels should be measured in a rest situation. This means that pumping from the well should be stopped a number of hours before taking the measurement. Yet in practice this is not always possible. Otherwise, when the stop is too short, the measured water level provides erroneous information. Hence, the recorded water levels are always examined and searched for outliers and other types of suspicious records. Once we find an erroneous data we delete it. Then it is necessary to substitute a ‘correct’ value instead of the wrong one.

The requirement for a complete set of data is typical to several applications. Consider, for example, the plotting of a contour map of the water levels surface. Recall that contour maps provide a useful tool for studying the behavior of an aquifer. In particular, a contour map enables us to draw stream lines, to compute gradients, and to detect the main flow directions. A comparison of two consecutive contour maps enables us to see how the groundwater surface changes from one year to another, and allows us to compute the relevant changes in the storage volume of the aquifer. Nowadays contour mapping is often done by a computerized plotting system whose input is a list of triplets, where each triplet contains the coordinates and the water level of one well. If one well is located far from the other wells then the water level in this well has a great effect on the shape of the computed water table surface. Deleting the well from the input list is likely, therefore, to cause a severe distortion to the shape of the computed surface. Similar distortions may occur when the missing water level is considerably higher (or lower) than the other water levels.

A further example is the case when the data are needed for the use of numerical models such as finite difference or finite element models. Here the observed water levels are used in the calibration process, to test the quality of the proposed model, and for setting initial and boundary conditions. In this case it is often necessary to transfer the data collected on a network of observation wells into a mesh of cells or elements which constitute the model. Usually each cell (or element) is associated with a small number of neighboring wells, and the water level at the cell is defined as a certain weighted average of the water levels in these wells. This averaging operation is often done by some computer program which does not take into account the effect of missing observations. So missing data are likely to cause either a breakdown or serious malformation to the output.

Missing water level observations are expected to cause similar difficulties in almost any monitoring system that is based on a permanent network of observation wells. The number of observation wells that constitute such a network can be quite large. For example, the Hydrological Service of Israel stores and updates records of over 6000 wells, which are used for monitoring the hydrological regime in several aquifers. In the past the majority of these wells were measured monthly. But twenty years ago the monitoring policy has been changed, and since then in most of the wells water levels are measured only twice a year: Once in the spring and once in the autumn. (However, a small percentage of important wells are still monitored on a monthly basis.) The reasons for this change of policy are two: first, a lack of manpower and budget; second, the belief that reports and models that are based on annual records of data are satisfactory for most purposes. (A similar policy is used in sampling concentrations of chloride and nitrate in wells. These samples are taken once a year, during the dry season.)

Note that there is substantial difference between the imputing of monthly records of water levels and the imputing of annual records. Assume for a moment that we have monthly records of water levels in a certain well, which were measured during m consecutive years, but contain a small percentage of missing records. Recall that monthly collected water levels are expected to show clear seasonal behavior, i.e. increasing during the rainy season and decreasing during the dry season. The seasonal pattern, and the fact that we have nearly m x 12 records per well, enable us to complete the missing records in a certain well by relying only on its own records (Dax 1985, and references therein). However, when having (at most) one record per year, the information gathered in one well is not sufficient for reliable imputing of the missing records. Consequently we are forced to rely on relevant information from neighboring wells, while the neighboring wells may also suffer from missing records. This raises the need for solving the above matrix completion problem.

The traditional treatment of missing data in hydrology is based on statistical methods. For example, autoregressive techniques for time series, or geostatistical techniques, like kriging, for spatial data. Nowadays statistical computing packages, like SAS, SPSS or R include built-in programs that achieve imputing. (For example, R users can choose between Amelia, Mice, and mitools.) Hence it is always tempting to use an off-the-shelf package. Indeed, in some cases we obtain reasonably good results. Yet, as explained in the next section, the statistical approach has certain limitations.

In this paper we propose a different approach. A ‘matrix imputing’ approach that takes advantage of matrix properties that can be derived from the matrix at hand. This approach has been proved very successful in several modern fields, such as machine learning, data mining, DNA microarrays data, computer vision, recommender systems, and collaborative filtering. However, as far as we know, it has not yet been tested in hydrology. One aim of this paper is, therefore, to introduce this approach, to explain the main ideas behind the method, and to illustrate how it works. This is achieved by describing the leading algorithms that characterize this method. Numerical experiments demonstrate the usefulness of this approach in the context of missing groundwater data.

The implementation of the above approach raises two issues. First we need an effective algorithm for computing a rank-k approximation. Second, we need to find a suitable value of k. The last problem is closely related to the question of how to determine the number of principle components (principle factors) of a data matrix, e.g., Howard & Gordon (1963), Jackson (1993), Peres-Neto et al. (2005) and Ledesma & Valero-Mora (2007). Yet, when A has missing entries we face a further difficulty: As k increases behind a certain threshold, the number of degrees of freedom in Problem (1) exceeds the number of known entries. At this point the solution of Problem (1) loses its uniqueness and the quality of the imputed values starts to deteriorate. A further aim of the paper is to illustrate this phenomenon and to outline an effective solution method.

The solution proposed in this paper is based on three innovations. The first one is an improved building process that generates a sequence of matrices, B₁, B₂, …, B_k, …, where B_k is a rank-k matrix that solves Problem (1). The building process is based on a new minimization procedure that uses B_k as a starting point for the computation of B_k+1. This results in a considerable saving of computing efforts. The second innovation regards the construction of a ‘virtual objective function’ that measures the quality of the imputed entries. This enables us to compute a sequence of positive numbers, V₁, V₂, …, V_k, …, where V_k reflects the ability of B_k to estimate the missing entries in A. The smaller is V_k, the better imputing we have. The third innovation is an effective rule for choosing an optimal value of k. As we shall see, the virtual objective function decreases in the first steps. Yet, as k increases behind a certain point, the values of V_k starts to increase. The optimal imputing is obtained, therefore, for the smallest value of V_k.

The plan of the paper is as follows. The Methodology section is divided into three parts. The first part explains the main differences between the traditional statistical approach and the current matrix approach. The second part outlines some of the leading matrix algorithms. There are several ways to compute a rank-k approximation. This can be done by direct solution of Problem (1), or by solving a closely related problem. In our review the focus is on methods that are most suitable for groundwater data. Starting from averaging methods we describe some basic imputing algorithms, including Iterative Column Regression (ICR), k – Nearest Neighbor (KNN) impute, and iterative Singular Value Decomposition (SVD) impute. Then we move to consider recently proposed methods such as rank minimization, and nuclear norm minimization. The third part of the Methodology section introduces a new matrix method that combines the three innovations that we mentioned above. The basic building block of the new method is a simple algorithm that computes a rank-one approximation of matrix. The usefulness of the proposed scheme is illustrated on a number of typical test cases.

The analysis of sample surveys is a major task in statistics. The problem of nonresponse in surveys is an inherent difficulty that hinders users from achieving accurate analysis. The statistical treatment of missing data is motivated, therefore, by the need to solve this problem, which leads to several differences between the statistical approach and the proposed matrix approach.

One feature that characterizes statistical imputing is that the methods rely on statistical assumptions on the origin of the data and the nature of the missing entries. For example, in maximum likelihood expectation-maximization (EM) methods, and in multiple imputing methods, it is common to assume that the entries of the matrix are generated by some known probabilistic distribution function (such as multivariate normal distribution). In contrast, matrix methods make no assumptions of this kind. Moreover, annual water levels are not expected to obey such assumptions.

Another type of assumption regards the missingness mechanism that determines the locations of the missing entries. Here it is common to distinguish between ‘not missing at random’, ‘missing at random’, and ‘missing completely at random’ (Dempster et al. 1977; Little & Rubin 1987; Allison 2009). Again, in matrix methods there is no need in such characterization, while missing water levels are not expected to obey these classifications.

The main difference between the statistical approach and the matrix approach lies in their goals.

Matrix methods have one aim: To provide ‘suitable’ values for the missing entries. In contrast, the ultimate aim of statistical analysis is to produce reliable estimates of statistical parameters that characterize the data, such as means, variances, correlations, covariances, principle components, and so on. If this task can be satisfied without assigning values to missing entries, then the analysis is carried out without achieving imputing. Consequently, ‘listwise deletion’ is the default option in most of the statistical software packages. In matrix terminology, listwise deletion means that we simply delete any row with missing entries. Clearly, this is not a reasonable solution for our problem. Yet the simplicity of this method makes it a popular option in many statistical applications. For example, as reported in Tsikriktsis (2005), an inspection of 103 research articles shows that listwise deletion is heavily used in operations management surveys. However, listwise deletion may lead to inaccurate estimates of statistical parameters, especially when the data fail to satisfy the ‘missing completely at random’ assumption. Thus, from the statistical point of view, the reason for conducting imputation is to gain improved estimates of the desired parameters. Moreover, in many methods the computed values are not guaranteed to help in other applications. For example, Allison (2009) gives the following warning about the use of imputed values that are obtained in the maximum likelihood method:

In statistical terminology, methods like mean substitution, hot-deck imputation, and EM, are classified as single imputation methods, since here each missing entry is assigned a single value. These methods have certain advantages over listwise deletion, but are liable to suffer from certain drawbacks. For example, mean substitution is likely to produce smaller estimates of the variance. The multiple-imputing approach is designed to avoid biased estimates (Little & Rubin 1987; Rubin 1987). In this method the values of the missing entries are randomly drawn from a given probabilistic distribution. Usually a multivariate normal distribution which is supplied by the user. Once all the missing entries in the matrix attain values, the algorithm computes the desired set of statistical parameters. This process of random draws and analysis is repeated ℓ times, generating ℓ matrices with ℓ sets of parameters. Typically ℓ is a small integer between 3 and 10. Then the final estimates of the parameters are obtained by taking means over the ℓ sets. This process generates an improved set of statistical parameters, as well as some indication on the accuracy of these parameters. Similarly, values for the missing entries can be obtained by averaging the corresponding entries in the ℓ matrices. However, as noted above, the assumption that water levels behave like random draws is somewhat doubtful.

Below we describe some popular algorithms and explain the main ideas behind these methods. The review is not extensive and many good methods have been left unmentioned. It is aimed to serve as a short introduction into the concept of matrix-imputing methods, with a focus on methods that are suitable for handling annual records of water levels.

All the methods share the same aim: To complete the missing entries of a given m x n data matrix, A = (a_ij). The description is based on standard terms and techniques of numerical computing. No assumptions are made on the statistical nature of the observed data.

Let α denote the average value of the known entries in A. The simplest way of imputing is, perhaps, substituting each missing entry of A by α. This mean-imputing method is useful in handling random matrices, whose entries are considered as random numbers that obey the same probabilistic distribution. However, the mean-imputing method is not suitable for groundwater data. Assume, for example, that a certain row of A refers to a rainy year. Then the (missing) water levels in this row are expected to exceed α. Similarly, consider the case when certain column of A refers to a well whose water levels are considerably lower than those of the other wells. Then the (missing) water levels in this column are expected to stay below the overall average.

Let y = (y₁, …, y_m)^T and w = (w₁, …, w_n)^T denote the corresponding vectors of unknowns. Let e = (1, 1, …, 1)^T denote a vector of ones whose length is either m or n, depending on the context. Recall that w denotes a column vector, w^T = (w₁, …, w_n) denotes a row vector, and the product y w^T defines an m x n rank-one matrix whose (i,j) entry equals y_i w_j, i = 1, …, m, j = 1, …, n. With these notations at hand the six averaging schemes mentioned above are essentially low-rank approximations of A regarding the matrices αe e^T, y e^T, e w^T, y e^T + e w^T, y e^T + e w^T +αe e^T, and y w^T, respectively. These observations suggest that higher rank approximations are likely to produce smaller model errors. Yet, averaging schemes are often used to compute starting values for the more sophisticated methods that are described below.

This method is most suitable to cases where the number of years, m, is considerably larger than the number of wells, n. As its name says, ICR is an iterative algorithm. It starts by setting initial trial values to the missing entries of A. Then these values are updated during the iterative process. The basic iteration is composed of n steps, where the jth step updates the values of the missing entries in the jth column, j = 1, …, n. This is achieved by regressing the jth column against the other columns of A. The details are given below.

The jth step of the basic iteration, j = 1, …, n.

• (a)

  Let denote the current m x n matrix at the beginning of the jth step.

• (b)

  Let the m x (n − 1) matrix A_j be obtained from by deleting the jth column of ⁠.

• (c)

  Let the m–vector b denote the jth column of ⁠.

• (d)
  Compute x*, the minimum norm solution of the linear least-squares problem

  
• (e)

  Compute the m-vector d = A_jx*.

• (f)

  Set new trial values to the missing entries in the jth column of A. The new values are those of the corresponding entries in d.

In statistician's eyes the above algorithm can be viewed as descendant of the statistical EM approach. Yet it is not clear whether the iterative process always converges, nor what properties a limit point has. Another related question is how the initial trial values of the unknown entries affect the limit point. The experiments in Hastie et al. (1999) reveal a slow rate of convergence (as is typical of the EM), but the imputed values are satisfactory. A similar behavior is observed in our experiments.

The ICR algorithm is based on the assumption that we have a small group of neighboring wells whose water levels behave in a similar way. So regressing one well against the others is likely to provide good results. However, if we have a large group of wells, from various places in the aquifer, there is no point in regressing one well against all the others. In this case the regression should be restricted to a smaller group of wells. This idea motivates the next method.

The KNN algorithm is a simple method that is able to handle a large number of wells. It relies on the assumption that for each well with missing entries we can find a small group of neighboring wells. The method is not iterative. It achieves only one sweep over the columns of A, considering one column at a time.

The treatment of the jth column, j = 1, …, n.

• (a)

  Let the m-vector b denote the jth column of A.

If b has no missing entries skip to the next column.

• (b)

  Select k other columns of A which are most ‘similar’ to b.

These columns are used to construct an m x k matrix, A_j.

• (c)

  Compute a k-vector d that estimates b. (Usually d is a weighted average of the k nearest neighbors.)

• (d)

  Set the missing entries in b to be the corresponding entries of d.

The above four steps are repeated for each column of A. Thus for each column we construct a different matrix of ‘neighboring’ columns, and a different simulating vector, d. Usually the number of neighbors, k, is determined in advance. A typical value of k lies between 4 and 10. It is also helpful to assign initial trial values to the missing entries. This simplifies the implementation of steps (b) and (c). For further discussions of this method see Hastie et al. (1999), Troyanskaya et al. (2001), and Kim et al. (2005).

The basic iteration:

• (a)

  Given an admissible matrix, A_ℓ, compute T_k(A_ℓ).

• (b)

  Compute a new admissible matrix A_ℓ+1. The entries of A_ℓ+1 which refer to missing entries in A obtain the values of the corresponding entries in T_k(A_ℓ).

The iterative SVD algorithm is another example of a method that descends from the EM framework. It alternates between the imputing process (expectation) and the SVD computation (maximization). The success of the method depends on proper choices of k and the starting matrix. A common way to compute A₁ is by applying an averaging scheme. It is shown in Dax (2014) that the related sequence of objective functions values is monotonic decreasing, which suggests that the sequence T_k(A_ℓ) may converge toward a solution of Equation (1.1). For further discussions of this method see Hastie et al. (1999), Troyanskaya et al. (2001), Friedland et al. (2006a, 2006b), and Kurucz et al. (2007).

Let u_ℓ-1 and v_ℓ-1 denote the current estimate of the solution at the beginning of the ℓ-th iteration. Then the ℓth iteration, ℓ = 1, 2, …, is composed of the following two steps.

Observe that minimizing f(u,v) by changing one variable at a time, results in the same basic iteration. Define f_ℓ = f (u_ℓ,v_ℓ) then, clearly, f_ℓ ≥ f_ℓ+1≥0 for ℓ = 1, 2, …, and the sequence {f_ℓ} converges.

The above algorithm is called Successive Rank-One Modifications (SROM). It is, perhaps, the simplest way to solve (1). Other methods for solving this problem include the Alternating Least Squares (ALS) algorithm, and the methods of Newton, Gauss-Newton, and Wiberg, e.g. Dax (2014, 2017) and Chen & Suter (2004). A major factor that effects the efficiency of the algorithm is the choice of the starting point. This issue is resolved below.

We have seen that F_k reflects the distance between known entries of A and the corresponding entries of B_k. It is tempting to think that F_k is also reflecting the ability of B_k to assess the unknown part of A. However, a small value of F_k does not guarantee that the entries of B_k provide a good estimate for the unknown entries of A. Moreover, as k increases behind a certain value, the number of degrees of freedom in (1) exceeds the number of known entries. At this point the solution of (1) loses its uniqueness, and the quality of the imputed values starts to deteriorate. (The number of degrees of freedom is close to k(m + n), which is the number of unknown variables in the presentation (16) of a rank-k matrix.) Yet, in spite of this deterioration, the values of F_k continue to decrease. This difficulty raises the need for another objective function that is able to assess the quality of the imputed entries.

However, in real life applications, the ‘true’ values of the missing entries are never known, so the definition of the virtual objective function needs some amendments. The modified definition is based on a small test set of entries which are chosen from the known entries of A. These entries are removed from the list of known entries and added to the list of unknown entries. Then the value of F_k in (17) is computed according to the new list. Similarly, V_k is computed by using (19) with the test set, while the value of μ in (20) gives the number of entries in the test set.

The idea of using formal and virtual objective functions is borrowed from the statistical Cross-Validation method. However, it seems to be a new idea in the context of groundwater imputing. The theory behind this method attempts to answer the questions of how many entries to include in the test set, and how to choose the test entries from among the known ones. In our experiments the entries of the test set are randomly chosen from the known entries of A.

Below we outline a simple rule for determining a value of k for which B_k achieves the best imputing. Once we decide on the value of k, the unknown entries of A are replaced by the corresponding entries of B_k. The basic tools for achieving this goal are the related sequences of formal and virtual objective functions values.

We have seen that the sequence of formal objective function values is always monotonic decreasing. In contrast, the virtual objective function values, V₁, V₂, …, stop to decrease after a few iterations. The ‘optimal’ value of k is that for which V_k attains its smallest value. Usually a brief inspection of the two sequences is sufficient for determining a suitable value of k.

In this section we present the results of some experiments with the proposed imputing methods. The algorithms were programmed in Fortran, using double precision arithmetic with round off unit about 10⁻¹⁶. The water level matrices are taken from the archive of the Hydrological Service of Israel and represent typical cases of missing records. As before, each data matrix, A = (a_ij), has m rows and n columns, where a_ij denotes the measured water level (if available) at the end of the ith hydrological year in the jth well. Thus each matrix contains the measured data during m years in a group of n neighboring wells. The names of the water level matrices come from the area in which the wells are located. The first three groups of wells, Rehovot, Shiflat Lod and Rishon, are located in the coastal aquifer of Israel. The other three groups, Ein Kerem, Yartan 4, and Yartan 5, are located in the Yarkon-Taninim (‘Yartan’) mountain aquifer of Israel.

The details of the data matrices are summarized in Table 2. The reading of this table is quite straightforward. Consider for example the third column, which considers the ‘Rishon’ matrix. Then this matrix has n = 6 columns, and m = 46 rows, while the number of missing entries, μ, equals 17. Note that the matrices Yartan 4 and Yartan 5 contain no missing entries.

The imputing methods that we have tested include three methods: the additive model, the ICR method, and the new algorithm, the Improved Building Process. The last method builds a sequence B₁, B₂, …, where B_k is a rank-k matrix that solves (1). (We have not tested the KNN imputing because it is mainly aimed at handling larger sets of wells.)

The data in Tables 3–5 provide the corresponding values of formal and virtual objective functions. The virtual objective function is defined by constructing a test set made of τ entries, which are randomly chosen from the known entries of A. The actual number of missing entries is, therefore, μ + τ. The other m_*n – (μ + τ) entries are used for constructing the formal objective function.

The matrix B₀ denotes the null matrix, a matrix with zero entries. The formal and the virtual values obtained for this matrix enable us to assess the improvement made by the tested methods. In additive model imputing the related formal and virtual values are computed by using the resulting rank-2 matrix, y e^T + e w^T + α e e^T, instead of B_k. The ICR method does not produce an approximating matrix and, therefore, only the related virtual values are computed. The starred figures denote minimal values of the virtual objective function. The smallest value of the formal objective function is always obtained for B_n.

The experiments with the first four test matrices are summarized in Table 3. Consider, for example, the experiments with the Rishon matrix. Running the additive model method on this matrix results in a formal objective value of 1.21, while the related value of the virtual objective function is 1.36. Similarly, when using the new method on the Rishon matrix, the matrix B₃ results in a formal objective function of 0.31 and a virtual objective function value of 0.63. The last value is the smallest virtual value obtained for the Rishon matrix. The ICR method gave a virtual value of 0.73, which is slightly above that value. The experiments in Table 3 clearly illustrate the usefulness of the proposed method for determining a value of k for which B_k achieves the best imputing.

The experiments in Tables 4 and 5 provide a further demonstration of the above advantage of the new method. In addition, these experiments enable us to see how the number of missing entries affects the imputing process. For this purpose we have used the matrices Yartan 4 and Yartan 5 with three different values of τ. (Recall that τ denotes the number of entries in the test set.) In these matrices μ = 0, so the overall number of missing entries equals τ. Let us consider, for example, the experiments on the Yartan 4 matrix when τ = 18. In this case the matrix B₂ results in a formal objective function value of 0.10 and a virtual objective function value of 0.23. The last value is a minimal virtual value, which is shared with the ICR algorithm. Similarly, let us consider the Yartan 5 matrix with τ = 16. In this case the minimal virtual value, 0.34, is obtained for B₂.

Graphical presentation of the imputed water levels are displayed in Figures 1 and 2. In these figures each curve refers to a different well, and the imputed values are marked with bold characters. The first figure considers the Ein-Kerem matrix, using B₃. (The measured water levels of the wells in Figure 1 are displayed in Table 1.) The second figure considers the Yartan 4 matrix with τ = 18, using B₂. These figures clearly illustrate the ability of the new method to achieve satisfactory imputing.

The current experiments use typical records of annual water levels, which were taken from the archive of the Hydrological Service of Israel. The statistical approach suggests further tests, with various scenarios of missing data. Such as ‘missing at random’ and ‘not missing at random’. However, this interesting issue is left for future research.

The need for imputing missing records of annual groundwater data arises in many applications. Traditional approaches for handling this problem are based on statistical techniques for treating missing data. In this paper we introduce a different approach, one that is based on matrix properties of the data matrix. The new approach has been found successful in many modern areas, but it has not yet been used in hydrology. Hence, one aim of this paper is to introduce the matrix approach and to demonstrate its usefulness. Indeed our experiments illustrate the ability of the new method to achieve satisfactory imputation of missing groundwater observations.

The method proposed in this paper is based on three important innovations. The first is a gradual rank increasing process that efficiently builds a sequence of matrices, B_k, k = 1, 2, 3, …, where B_k is a rank-k matrix that solves (1). The second innovation regards the construction of a ‘virtual objective function’ that measures the quality of the imputed entries. This enables us to build a related sequence, V_k, k = 1, 2, 3, …, where V_k reflects the ability of B_k to estimate the missing entries in A. The third innovation is an effective rule for choosing an optimal value of k, for which B_k achieves the best imputing. As we have seen, the virtual values, V₁, V₂, …, stop decreasing after a few iterations, and the optimal value of k is that for which V_k attains its smallest value. Another benefit of using a ‘virtual objective function’ is that it enables us to compare the performance of different methods. This point is illustrated in our experiments, which compare the additive model, the ICR algorithm, and the new method.

The authors are greatly indebted to three anonymous referees for several helpful comments.