Some Algorithms for the Conditional Mean Vector and Covariance Matrix

We consider here the problem of computing the mean vector and covariance matrix for a conditional normal distribution, considering especially a sequence of problems where the conditioning variables are changing. The sweep operator provides one simple general approach that is easy to implement and update. A second, more goal-oriented general method avoids explicit computation of the vector and matrix, while enabling easy eval-uation of the conditional density for likelihood computation or easy generation from the conditional distribution. The covariance structure that arises from the special case of an ARMA( p, q ) time series can be exploited for substantial improvements in computational eﬃciency.

Note that the conditional density of (Z 1 | Z 2 = z 2 ) equals the ratio of the joint normal density of (Z 1 , Z 2 ) and the marginal normal density of Z 2 evaluated at Z 2 = z 2 . Taking µ = 0 to simplify, cancelling constants, and taking the log of both sides, we get the useful result which can be rewritten as Amidst the cancellation of the constants is a useful relationship of the determinants: We consider here the problem of computing the conditional mean vector µ 1.2 and covariance matrix A 1.2 where n is large, s is modest in size, and where the partitioning is arbitrary. The applications envisioned include both generation from the conditional distribution and inference in cases of censored or missing values. Especially interesting are applications where the partitioning may be repeatedly changing, as may occur in some MCMC techniques. In Sections 2 and 3, we consider the case of the general multivariate normal distribution, and the special case where Z arises from a Gaussian autoregressive-moving average (ARMA(p, q)) process is considered in Section 4. FORTRAN 95 code implementing some of these results is included with this manuscript as 'gsubex1.f95' to 'gsubex3.f95'.

A general approach with the sweep operator
Following standard computational techniques (see, for example, Golub and van Loan 1996;Monahan 2001;Stewart 1973) such as Gaussian elimination or Cholesky factorization, the direct computation of µ 1.2 and A 1.2 will take O(n 3 ) floating point operations, owing from computing the inverse (solving a system of equations) in A 22 which is O(n) in size. However, if the partitioning were to be changed slightly, little of these results could be used to address the modified problem. The sweep operator is more appropriate here, being easy to update while taking nearly the same computations to form µ 1.2 and A 1.2 .
The sweep operator is widely used in regression problems (see Stiefel 1963;Goodnight 1979;Monahan 2001, Section 5.12) computing the following tableau after sweeping the first p rows/columns of the (p + 1) × (p + 1) matrix: The sweep operator can be viewed as a compact scheme for full column elimination using diagonal pivots v jj . It is usually implemented as a single step, operating on one row/column at a time, so that the inverse is a reciprocal. When an identity matrix is augmented, as below, full elimination pivoting on the first p columns produces the matrix Merely overwriting the inverse on top of the matrix, and avoiding storing the known elements of the identity matrix produces the sweep step. The sweep operator can be used to compute the regression coefficients (X T X) −1 X T y and error sum of squares y T y−y T X(X T X) −1 X T y. Applying the sweep operator to the partitioned covariance matrix A, the result of sweeping the LAST n − s rows/columns produced the following tableau which obviously includes the necessary information for the conditional distribution. The ease of computing updates follows from the commutativity and reversibility properties of the sweep operator. The commutativity property means that the order of the sweeping does not matter. Reversibility means that sweeping a particular row/column a second time is the same as not sweeping it in the first place. As a consequence, sweeping a particular row/column merely moves a variable from one subset to the other. For the conditional distribution, the rows/columns corresponding to the variables that are to be conditioned upon are swept.
Conditioning on a new variable means sweeping that row/column of the covariance matrix. Removing a variable from the list to be conditioned upon means sweeping that row/column a second time. Each sweep of a row/column takes O(n 2 ) operations.
As in computing the determinant from full column pivoting, the determinant of A 22 can be computed from the product of the pivot elements v jj , which can be seen in the following Example 1.
Example 1: Let n = 3, and A given below and the following sweep sequence: Note that the repeating decimals .333 and .166 in the example above are rounded versions of 1/3 and 1/6. In practice, rounding error will accumulate with each sweep step, and the process may need to be restarted when the accumulated error becomes noticeable.

A second general method
This second general approach aims to solve the two specific problems without explicitly computing the conditional mean vector and covariance matrix. One task is computing the likelihood of the conditional normal distribution, namely the quadratic form and the determinant |A 1.2 |. The second task would be to generate a random vector from the conditional distribution.
Consider the reverse-Cholesky factorization of the covariance matrix A, producing the upper triangular factor U : The usual Cholesky factorization produces a lower triangular matrix and the induction order goes from upper left to lower right; here the order is reversed and an upper triangular matrix is produced. Now consider the problem of generating from the joint distribution using a random vector v ∼ Normal n (0, I n ).
Partitioning v in the same way, the algebra would follow the simple form Conditioning on Z 2 = z 2 would really mean knowing z 2 . If that were so, we could then solve for v 2 , solving the triangular system to get v 2 = U −1 3 (z 2 − µ 2 ).
Then z 1 could then be generated conditional on the value of z 2 , following leading to the result A 1.2 = U 1 U T 1 . For computing the likelihood, the determinant part is easy: |A 1.2 | = | U 1 | 2 , following from (5), since | U 3 | 2 = |A 22 |. The quadratic form is nearly as easy, since from (4) we have The quadratic form in A −1 can be written in partitioned form as the squared length of , while the quadratic form in (A 22 ) −1 is just the squared length of the second partitioned vector above so that the difference is the squared length of the first component As presented, this approach offers no particular advantage over any other direct method -it merely uses a form of Cholesky decomposition. If the partitioning were changed, however, this Cholesky factor can be updated. The updating procedure due to Clarke (1981) was devised for subset regression computations and follows the same spirit as the sweep operator (see also Daniel, Gragg, Kaufman, and Stewart 1976). The efficiency of the update scheme will depend on the pattern of changes of partitioning.
The change in partitioning can be viewed in terms of a permutation matrix P . Conditioning on an additional component (or one fewer) means permuting that component to the partition boundary and then moving the boundary. Here we want to construct the reverse-Cholesky (upper triangular) factor of the changed matrix P AP T using P U . The problem is that P U is no longer upper triangular but needs to be returned to that form. This can be done by postmultiplying by sufficient Givens transformations G 1 , G 2 , ..., G m , each an orthogonal matrix, as many as may be needed to return P U G 1 · · · G m to upper triangular form, and restoring the factorization Example 2: For the covariance matrix A below, the decomposition provides U that could be easily employed for (Z 1 , Z 2 |Z 3 = z 3 , Z 4 = z 4 ). Now suppose we now To return this matrix to upper triangular form, first postmultiply by a Givens transformation G 1 defined by elements (4,3) and (4,4), changing everything in columns 3 and 4, and placing a zero in (4,3). Here G 1 is an identity matrix except for (G 1 ) 33 = (G 1 ) 44 = 1 √ 5 and (G 1 ) 34 = −(G 1 ) 43 = 2 √ 5. Next construct G 2 using (3,2) and (3,3), changing columns 2 and 3 and placing a zero in (3,2), and return the upper triangular form. The zeros in (2,2) and (2,3) will change with these two transformations, but this will not affect the triangular form.

ARMA covariance model
The case where Z t , t = 1, . . . , n follows the ARMA(p, q) time series model brings some special structure that can be exploited. If the partitioning followed the same simple structure as in (1), then the conditional distribution can be easily constructed. However, we are interested here in the case where the partitioned vectors are Z 1 = (Z S(1) , Z S(2) , . . . , Z S(s) ) T and Z 2 = (Z T (1) , Z T (2) , . . . , Z T (n−s) ) T and where subset list of indices S(j) may be arbitrary and S and T are complementary.
A device due to Ansley can be exploited to compute the conditional distribution efficiently. Let Z t , t = 1, . . . , n follow the ARMA(p, q) process (6) and denote the covariance matrix by A. Now let the n × n matrix B have its first p rows all zeros except for 1 on the diagonal, and the last n − p rows of the form aligned so that the 1 appears on the diagonal. Ansley (1979) shows that the covariance matrix of BZ, namely BAB T , is a banded positive definite matrix. The bandwidth of BAB T is max(p, q) = m for the first p rows, and q thereafter. Its Cholesky factor L, such that is lower triangular with the same banding structure, and can be computed in O(nm 2 ). Solving a system of equations in L, say to compute L −1 u, can done in O(nm). The space required is only O(m 2 ). As a result, a bilinear form in the inverse of the covariance matrix can be computed in O(nm 2 ) time and O(m 2 ) space.
Before applying this result to the arbitrary partitioning problem, and so following the original partitioning in (1), first notice that which can be established through the usual results for the inverse of a partitioned matrix. The remaining results are more complicated. Defining Q(z) = z T A −1 z, and employing result (4), we have the derivative results ∂ 2 Q(z)/∂z i ∂z j = 2 (A 1.2 ) −1 ij for 1 ≤ i, j ≤ s, and for 1 ≤ i ≤ s. For computation, we replace the derivative with a difference expression. In the case of a quadratic function Q(z), the second difference will be exact for the second partial derivative, so that where e i is the i th elementary vector, verifying (7) above. A similar result using the first difference, at z 1 = 0, gives Now to transform these results to the case of arbitrary parititioning, construct the n × s matrix P , which is zero everywhere except for elements (P ) S(i),i = 1 for 1 ≤ i ≤ s, and, similarly, the n × (n − s) matrix Q, which is zero everywhere except (Q) T (i),i = 1 for 1 ≤ i ≤ n−s. Putting P and Q together as P Q forms an n×n permutation matrix. Then we have the computational formulas The key computational advantage of this method is that, as bilinear forms in A −1 , the computational burden is only O(nm 2 + nms + s 2 ) for ARMA(p, q) time series models.
Finally, if the conditioning variables z 2 (or its mean) were changed, the resulting vector u would have to be recomputed, taking only O(n). If a variable were moved from S to T , the mean vector would change and u must be updated. If a variable were moved from T to S, then a new column would be added to P , and the mean vector would also have to be updated.