{'fraction_non_alphanumeric': 0.04588185772414761, 'fraction_numerical': 0.013750526168093167, 'mean_word_length': 4.051204819277109, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Block BFGS: Squeezing More Curvature out of Data

We propose a novel limited-memory stochastic block BFGS update for incorporating enriched curvature information in stochastic approximation methods. In our method, the estimate of the inverse Hessian matrix that is maintained by it, is updated at each iteration using a sketch of the Hessian, i.e., a randomly generated compressed form of the Hessian. We propose several sketching strategies, present a new quasi-Newton method that uses stochastic block BFGS updates combined with the variance reduction approach SVRG to compute batch stochastic gradients, and prove linear convergence of the resulting method. Numerical tests on large-scale logistic regression problems reveal that our method is more robust and substantially outperforms current state-of-the-art methods.


Introduction
We design a new stochastic variable-metric (quasi-Newton) method-the stochastic block BFGS method-for solving the Empirical Risk Minimization (ERM) problem: We assume the loss functions f i : R d → R to be convex and twice differentiable and focus on the setting where the number of data points (examples) (n) is very large.
To solve (1), we employ iterative methods of the form where η > 0 is a stepsize, g t ∈ R d is an estimate of the gradient ∇f (x t ) and H t ∈ R d×d is a positive definite estimate of the inverse Hessian matrix, that is H t ≈ ∇ 2 f (x t ) −1 . We refer to H t as the metric matrix 1 .
The most successful classical optimization methods fit the format (2), such as gradient descent (H t = I), Newton's method H t = ∇ 2 f (x t ) −1 , and the quasi-Newton methods H t ≈ ∇ 2 f (x t ) −1 ; all with g t = ∇f (x t ). The difficulty in our setting is that the large number of data points makes the computational costs of a single iteration of these classical methods prohibitively expensive.
To amortize these costs, the current state-of-the-art methods use subsampling, where g t and H t are calculated using only derivatives of the subsampled function where S ⊆ [n] def = {1, 2, . . . , n} is a subset of examples selected uniformly at random. Using the subsampled gradient ∇f S (x) as a proxy for the gradient is the basis for (minibatch) stochastic gradient descent (SGD), but also for many successful variance-reduced methods [32,36,18,19,5,35,20] that make use of the subsampled gradient in calculating g t .
Recently, there has been an effort to calculate H t using subsampled Hessian matrices ∇ 2 f T (x t ), where T ⊆ [n] is sampled uniformly at random and independently of S [6,31]. The major difficulty in this approach is that calculating ∇ 2 f T (x t ) can be computationally expensive and, when d is large, storing ∇ 2 f T (x t ) in memory can be infeasible. One fairly successful solution to these issues [2,25] is to use a single Hessian-vector product ∇ 2 f T (x t )v, where v ∈ R d is a suitably selected vector, to update H t by the limited-memory (L-BFGS) [26] version of the classical BFGS [1,7,9,37] method. Calculating this Hessian-vector product can be done inexpensively using the directional derivative In particular, when using automatic differentiation techniques [4,16] or backpropagation on a neural network [27], evaluating the above costs at most five times as much as the cost of evaluating the subsampled gradient ∇f T (x t ).
Using only a single Hessian-vector product to update H t yields only a very limited amount of curvature information, and thus may result in an ineffective metric matrix. The block BFGS method addresses this issue. The starting point in the development of the block BFGS method is the simple observation that, ideally, we would like the metric matrix H t to satisfy the inverse equation H t ∇ 2 f T (x t ) = I, since then H t would be the inverse of an unbiased estimate of the Hessian. But solving the inverse equation is computationally expensive. So instead, we propose that H t should instead satisfy a sketched version of this equation, namely where D t ∈ R d×q is a randomly generated matrix which has relatively few columns (q d). The sketched subsampled Hessian ∇ 2 f T (x t )D t can be calculated efficiently through q directional derivatives of the form (3).
Note that (4) has possibly an infinite number of solutions, including the inverse of ∇ 2 f T (x t ). To determine H t uniquely, we maintain a previous estimate H t−1 ∈ R d×d and project H t−1 onto the space of symmetric matrices that satisfy (4). The resulting update, applied to H t−1 to arrive at H t , is the block BFGS update.
In the remainder of the paper we describe the block BFGS update, propose a new limitedmemory block BFGS update and introduce several new sketching strategies. We conclude by presenting the results of numerical tests of a method that combines the limited-memory block BFGS update with the SVRG/S2GD method [18,19], and demonstrate that our new method yields dramatically better results when compared to SVRG, or the SVRG method coupled with the classical L-BFGS update as proposed in [25].

Contributions
This paper makes five main contributions: (a) New metric learning framework. We develop a stochastic block BFGS update 2 for approximately tracking the inverse of the Hessian of f . This technique is novel in two ways: the update is more flexible than the traditional BFGS update as it works by employing the actions of a subsampled Hessian on a set of random vectors rather than just on a single deterministic vector, as is the case with the standard BFGS. That is, we use block sketches of the subsampled Hessian.
(b) Stochastic block BFGS method. Our block BFGS update is capable of incorporating enriched second-order information into gradient based stochastic approximation methods for solving problem (1). In this paper we illustrate the power of this strategy in conjunction with the strategy employed in the SVRG method of [18] for computing variance-reduced stochastic gradients. We prove that the resulting combined method is linearly convergent and empirically demonstrate its ability to substantially outperform current state-of-the-art methods.
(c) Limited-memory method. To make the stochastic block BFGS method applicable to largescale problems, we devise a new limited-memory variant of it. As is the case for L-BFGS [26], our limited-memory approach allows for a user-defined amount of memory to be set aside. But unlike L-BFGS, our limited-memory approach allows one to use the available memory to store more recent curvature information, encoded in sketches of previous Hessian matrices. This development of a new limited block BFGS method should also be of general interest to the optimization community.
(d) Factored form. We develop a limited-memory factored form of the block BFGS update. While factorized versions of standard quasi-Newton methods were proposed in the 1970's (e.g., see [8,10]), as far as we know, no efficient limited memory versions of such methods have been developed. Factored forms are important as they can be used to enforce positive definiteness of the metric even in the presence numerical imprecision. Furthermore, we use the factored form in calculating new sketching matrices.
(e) Adaptive sketching. Not only can sketching be used to tackle the large dimensions of the Hessian, but it can also simultaneously precondition the inverse equation (4). We present a selfconditioning (i.e., adaptive) sketching that makes use of the efficient factored form of the block BFGS method developed earlier. We also present a sketching approach based on using previous search directions. Our numerical tests show that adaptive sketching can in practice lead to a significant speedup in comparison with sketching from a fixed distribution.

Background and Related Work
The first stochastic variable-metric method developed that makes use of subsampling was the online L-BFGS method [34]. In this work the authors adapt the L-BFGS method to make use of subsampled gradients, among other empirically verified improvements. The regularized BFGS method [24,23] also makes use of stochastic gradients, and further modifies the BFGS update by adding a regularizer to the metric matrix.
The first method to use subsampled Hessian-vector products in the BFGS update, as opposed to using differences of stochastic gradients, was the SQN method [2]. Recently, [25] propose combining SQN with SVRG. The resulting method performs very well in numerical tests. In our work we combine a novel stochastic block BFGS update with SVRG, and prove linear convergence. The resulting method is more versatile and superior in practice to SQN as it can capture more useful curvature information.
The update formula, that we refer to as the block BFGS update, has a rather interesting background. The formula first appeared in 1983 in unpublished work by Schnabel [33] on designing quasi-Newton methods that make use of multiple secant equations. Schnabel's method requires several modifications that stem from the lack of symmetry and positive definiteness of the resulting update. Later, and completely independently, the block BFGS update appears in the domain decomposition literature [22] as a preconditioner, where it is referred to as the balancing preconditioner. In that work, the motivation and derivation are very different from those used in the quasi-Newton literature; for instance, no variational interpretation is given for the method. The balancing preconditioner was subsequently taken out of the PDE context and tested as a general purpose preconditioner for solving a single linear system and systems with a changing right hand side [15]. Furthermore, in [15] the authors present a factored form of the update in a different context, which we adapt for limited-memory implementation. Finally, and again independently, a family of block quasi-Newton methods that includes the block BFGS is presented in [13] through a variational formulation and in [17] using Bayesian inference.

Stochastic Block BFGS Update
The stochastic block BFGS update, applied to H t−1 to obtain H t , is defined by the projection where and Tr (·) denotes the trace. The method is stochastic since D t ∈ R d×q is a random matrix. The constraint in (5) serves as a fidelity term, enforcing that H t be a symmetric matrix that satisfies a sketch of the inverse equation (4). The objective in (5) acts as a regularizer, and ensures that the difference between H t and H t−1 is a low rank update. The choice of the objective is special yet in another sense: recent results show that if the iteration (5) is applied to a fixed invertible matrix A in place of the subsampled Hessian, then the matrices H t converges to the inverse A −1 at a linear rate [13]. This is yet one more reason to expect that in our setting the matrices H t approximately track the inverse Hessian. For similar techniques applied to the problem of solving linear systems, we refer to [12,14]. The solution to (5) is This solution was given in [11,17] and in [33] for multiple secant equations.
Note that stochastic block BFGS (6) yields the same matrix H t if D t is replaced by any other matrixD t ∈ R d×q such that span(D t ) = span(D t ). It should also be pointed out that the matrix H t produced by (6) is not what would be generated by a sequence of q rank-two BFGS updates based on the q columns of D t unless the columns of D t are ∇ 2 f T (x t )-conjugate. Note that these columns would be conjugate if they were generated by the BFGS method applied to the problem of minimizing a strictly convex quadratic function with Hessian ∇ 2 f T (x t ), using exact line-search.
We take this opportunity to point out that stochastic block BFGS can generate the metric used in the Stochastic Dual Newton Ascent (SDNA) method [29] as a special case. Indeed, when H t−1 = 0, then H t is given by When the sketching matrix D t is a random column submatrix of the identity, then (7) is the positive semidefinite matrix used in calculating the iterates of the SDNA method [29]. However, SDNA operates in the dual of (1).

Stochastic Block BFGS Method
The goal of this paper is to design a method that uses a low-variance estimate of the gradient, but also gradually incorporates enriched curvature information. To this end, we propose to combine the stochastic variance reduced gradient (SVRG) approach [18] with our novel stochastic block BFGS update, described in the previous section. The resulting method is Algorithm 1. Algorithm 1 has an outer loop in k and an inner loop in t. In the outer loop, the outer iterate w k ∈ R d and the full gradient ∇f (w k ) are computed. In the inner loop, both the estimate of the gradient g t and our metric H t are updated using the SVRG update and the block BFGS update, respectively.
To form the sketching matrix D t we employ one of the following three strategies: (a) Gaussian sketch. D t has standard Gaussian entries sampled i.i.d at each iteration. (c) Self-conditioning. We sample C t ⊆ [d] uniformly at random and set D t = L t−1 I :Ct = [L t−1 ] :Ct , where L t−1 L T t−1 = H t−1 and I :Ct denotes the concatenation of the columns of the identity matrix indexed by a set C t ⊂ [d]. Thus the sketching matrix is formed with a random subset of columns of a factored form of H t . The idea behind this strategy is that the ideal sketching matrix should be D t = (∇ 2 f T (x t )) −1/2 so that the sketch not only compresses but also acts as a preconditioner on the inverse equation (4). It was also shown in [13] that this choice of sketching matrix can accelerate the convergence of H t to the inverse of a fixed matrix. In Section 3.2 we describe in detail how to efficiently maintain and update the factored form.

Limited-Memory Block BFGS
When d is large, we cannot store the d × d matrix H t . Instead, we store M block triples, consisting of previous block curvature pairs and the inverse of their products With these triples we can form the H t operator implicitly by using a block limited-memory two loop recurrence. To describe this two loop recurrence, Since we do not store H t for any t, we do not have access to H t−M . In our experiments we simply set H t−M = I (the identity matrix). Other, more sophisticated, choices are possible, but we do not explore them further here. Using the above expansion, the action of the operator H t on a vector v can be efficiently calculated using Algorithm 2.
Algorithm 2 Block L-BFGS Update (Two-loop Recursion) The total cost in floating point operations of executing Algorithm 2 is M q(4d + 2q). In our experiments M = 5 and q will be orders of magnitude less than d, typically q ≤ √ d. Thus the cost of applying Algorithm 2 is approximately O(d 3/2 ). This does not include the cost of computing the product D t T Y t (O(q 2 d) operations) and its Cholesky factorization (O(q 3 ) operations), which is done outside of Algorithm 2. The two places in Algorithm 2 where multiplication by ∆ i is indicated is in practice performed by solving two triangular systems using the Cholesky factor of D i T Y i . We do this because it is more numerically stable than explicitly calculating the inverse matrix ∆ i .

Factored Form
Here we develop a new efficient method for maintaining and updating a factored form of the metric matrix. This facilitates the development of a novel idea which we call self-conditioning sketch.
t . An update formula for the factored form of H t , i.e., for L t for which H t = L t L T t , was recently given (in a different context) in [15]: This factored form of H t is too costly to compute because it requires inverting L t−1 . However, if we let D t = L t−1 I :Ct , where C t ⊂ [d], then (9) reduces to which can be computed efficiently. Furthermore, this update of the factored form (10) is amenable to recursion and can thus be expanded as By storing M previous curvature pairs (8) and additionally the sets C t+1−M , . . . , C t , we can calculate the action of L t on a matrix V ∈ R d×q by using (11), see Algorithm 3. To the best of our knowledge, this is the first limited-memory factored form in the literature. Since we do not store L t explicitly for any t we do not have access to L t−M , required in computing (11). Thus we simply use L t−M = I.
Again, we can implement a more numerically stable version of Algorithm 3 by storing the Cholesky factor of D T i Y i and using triangular solves, as opposed to calculating the inverse matrix

Convergence
In this section we prove that Algorithm 1 converges linearly. Our analysis relies on the following assumption, and is a combination of novel insights and techniques from [19] and [25].
for all x ∈ R d and all T ⊆ [n].
We need two technical lemmas, whose proofs are given in Sections 7 and 8 at the end of the paper.
A proof of this lemma is given in Appendix 2, where in particular it is shown that the lower bound satisfies γ ≥ 1 1+M Λ and the upper bound satisfies where κ def = Λ/λ. We now state a bound on the norm of the SVRG variance-reduced gradient for minibatches. Lemma 2. Suppose Assumption 1 holds, let w * be the unique minimizer of f and let w, x ∈ R d . Let µ = ∇f (w) and g = ∇f S (x) − ∇f S (w) + µ. Taking expectation with respect to S, we have The following theorem guarantees the linear convergence of Algorithm 1.
Theorem 1. Suppose that Assumption 1 holds. Let w * be the unique minimizer of f . When Option II is used in Algorithm 1, we have for all k ≥ 0 that where the convergence rate is given by assuming we have chosen η < γλ/(2Γ 2 Λ 2 ) and that we choose m large enough to satisfy 3 m ≥ 1 2η (γλ − ηΓ 2 Λ(2Λ − λ)) .
Proof. Since g t = ∇f St (x t ) − ∇f St (w k ) + µ and x t+1 = x t − ηH t g t in Algorithm 1, from (12) we have that .
Taking expectation conditioned on x t (i.e., with respect to S t , T t and D t ) and using Lemma 1 we have Introducing the notation δ f (x) def = f (x) − f (w * ) and applying Lemma 2 and the fact that strongly convex functions satisfy the inequality ∇f (x) 2 2 ≥ 2λδ f (x) for all x ∈ R d , gives where α = 2η γλ − ηΓ 2 Λ 2 and β = 2η 2 Γ 2 Λ(Λ − λ). Taking expectation, summing over t = 0, . . . , m − 1 and using telescopic cancellation gives 3 By our assumption on η, the expression on the right is nonnegative.

9
where we used that w k = x 0 and m t=1 E [x t ] = mE [w k+1 ] which is a consequence of using Option II in Algorithm 1. Rearranging the above gives where we used that f (w * ) ≤ f (x m ). Using that η < γλ/(2Γ 2 Λ 2 ), it follows that

Numerical Experiments
To validate our approach, we compared our algorithm to SVRG [18] and the variable-metric adaption of SVRG presented in [25], which we refer to as the MNJ method. We tested the methods on seven empirical risk minimization problems with a logistic loss and L2 regularizer, that is, we solved where A = [a 1 , . . . , a n ] ∈ R d×n and y ∈ {0, 1} n are the given data. We have also employed the standard "bias trick" 4 . For our experiments we used data from the LIBSVM collection [3]. All the methods were implemented in MATLAB. All the code for the experiments can be downloaded from http://www.maths.ed.ac.uk/~prichtar/i_software.html. We tested three variants of Algorithm 1, each specified by the use of a different sketching matrix. In Table 1 we present a key to the abbreviations used in all our figures. The first three methods, gauss, prev and fact, are implementations of the three variants (a), (b) and (c), respectively, of Algorithm 1 using the three different sketching methods discussed at the start of Section 3. In these three methods the M stands for the number of stored curvature pairs (8) used.
For each method and a given parameter choice, we tried the stepsizes η ∈ {10 0 , 5 · 10 −1 , 10 −1 , . . . , 10 −7 , 5 · 10 −8 , 10 −8 } and reported the one that gave the best results. We used the error f (x t ) − f (w * ) for the y-axis in all our figures 5 . We calculated f (w * ) by running all methods for 30 passes over the data, and then taking the minimum function value. Finally, we used m = n/|S t | for the number of inner iterations in all variants of Algorithm 1, SVRG and MNJ, so that all methods perform an entire pass over the data before recalculating the gradient.

Parameter investigation
In our first set of tests we explored the parameter space of the prev variant of Algorithm 1. By fixing all parameters but one, we can see how sensitive the prev method is to the singled-out parameter, but also, build some intuition as to what value or, interval of values, yields the best results for this parameter. We focused our tests on the prev method since it proved to be overall, the most robust method.
In Figure 1 (a) we depict the results of varying the memory parameter M , while fixing the remaining parameters. In particular, we fixed |T t | = |S t | = 15. In both the subplots of error × datapasses and error × time in Figure 1 (a), we see that 1 ≤ M ≤ 4 resulted in the best performance. Furthermore, the error × datapasses is insensitive to increasing the memory, since increasing the memory parameter does not incur in any additional data passes. On larger dimensional problems we found that approximately M = 5 yielded the overall best performance. Thus we used M = 5 in our tests on large-scale problems in Sections 5.2 and 5.3.
In Figure 1 (b) we experimented with varying |T t |, the size of the Hessian subsampling. Here the results of the error × datapasses subplot conflict with those of the the error × time subplots. While in the error × time subplot the method improves as |T t | increases, in the error × datapasses subplot the method improves as |T t | decreases. As a compromise, in Sections 5.2 and 5.3 we use |S t | = |T t | as our default choice.
In Figure 2 (a) we experimented with varying the size of gradient subsampling and Hessian subsampling jointly with |S t | = |T t |. In both the error × datapasses and the error × time subplot the range n 1/3 = 32 ≤ |S t | = |T t | ≤ 724 = n 1/2 , resulted in a good performance. Based on this experiment, we use |S t | = |T t | = n 1/2 as our default choice in Sections 5.2 and 5.3. Note that when |S t | is large the method passes through the data in fewer iterations and consequentially less time. This is why the method terminated early (in time) when |S t | is large. Finally in Figure 2 (b) we vary the parameter L, that is, the number of previous search directions used to form the columns of D t . From these tests we can see that using a value that is too small, e.g. L = 1, or a value that is too large L = 2 √ d = 24, results in an inefficient method in terms of both time and datapasses. Instead, we get the best performance when d 1/4 ≤ L ≤ d 1/2 . Thus on the first five problems we used either L = d 1/3 or L = d 1/4 , depending on which gave the best performance. As for the last two problems: rcv1-train.binary and url-combined, we found that L = d 1/4 was too large. Instead, we probed the set L ∈ {2, . . . , 10} for an L that resulted in a reasonable performance.
Through these four experiments in Figures 1 and 2, we can conclude that the prev method is not overly sensitive to the choice of these parameters. That is, the method works well for a range of parameter choices. This is in contrast with choosing the stepsize parameter, whose choice can make the difference between a divergent method and a fast method.

Data passes
We now compare all the methods in Table 1 in terms of error × datapasses. Since our experiment in Figure 1 (b) indicated that |S t | = |T t | resulted in a reasonable method, for simplicity, we used the same subsampling for the gradient and Hessian in all of our methods; that is S t = T t . This is not necessarily an optimal choice. Furthermore, we set the subsampling size |S t | = √ n. For the MNJ method, we used the suggestion of the authors of both [2] and [25], and chose L = 10 and |T t | ≈ L|S t | so that the computational workload performed by L inner iterations of the SVRG method was approximately equal to that of applying the L-BFGS metric once. The exact rule we found to be efficient was We set the memory to 10 for the MNJ method in all tests, which is a standard choice for the L-BFGS method.
On the problems with d significantly smaller then n, such as Figures 1(c) and 1(d), all the methods that make use of curvature information performed similarly and significantly better than the SVRG method. The prev method proved to be the best overall and the most robust, performing comparably well on problems with d n, such as Figures 1(c) and 1(d), but was also the most efficient method in Figures 1(a), 1(b), 1(e) (shared being most efficient with MNJ) and Figure 2(a). The only problem on which the prev method was not the most efficient method was on the url-combined problem in Figure 2(b), where the MNJ method proved to be the most efficient.
Overall, these experiments illustrate that incorporating curvature information results in a fast and robust method. Moreover, the added flexibility of the block BFGS update to incorporate more curvature information, as compared to a single Hessian-vector product in the MNJ method, can significantly improve the convergence of the method, as can be seen in Figures 1(a) and 1(b).

Timed
In Figures 5 and 6, we compare the evolution of the error over time for each method. While measuring time is implementation and machine dependent, we include these time plots as to provide further insight into the methods performance. Note that we did not use any sophisticated implementation tricks such as "lazy" gradient updates [19], but instead implemented each method as originally designed so that the methods can be compared on a equal footing.
The results in these tests corroborate with our conclusions in Section 5.

Extensions
This work opens up many new research avenues. For instance, sketching techniques are increasingly successful tools in large scale machine learning, numerical linear algebra, optimization and computer science. Therefore, one could employ a number of new sketching methods in the block BFGS method, such as the Walsh-Hadamard matrix [28,21]. Using new sophisticated sketching methods, combined with the block BFGS update, could result in even more efficient and accurate estimates of the underlying curvature. Viewed in terms of sketching, the MNJ method [25] and our Stochastic Block BFGS method with prev sketching follow opposite strategies. While the MNJ method sketches the Hessian matrix with a single vector v, where v is the average over a combination of previous search directions (a very coarse approximation), our prev variant uses all previous search directions to form the sketching matrix D t (a finer approximation). Deciding between these two extremes or how to combine them could be done adaptively by examining the curvature matrix D t ∇f (x t )D t . When the previous search directions are almost collinear, D t ∇f (x t )D t becomes ill-conditioned. In this case one should form matrix D t with less columns using fewer or a coarser combination of previous search directions, while if D t ∇f (x t )D t is well-conditioned, one could use more or a finer combination of previous search directions.
While in this work we have for simplicity focused on utilizing our metric learning techniques in conjunction with SVRG, they can be used with other optimization algorithms as well, including SGD [30], SDCA [36] and more.

Proof of Lemma 2
Let h S (w) = f S (w) − f S (w * ) − ∇f S (w * ) T (w − w * ). Note that h S (w) achieves its minimum at w * and h S (w * ) = 0. Furthermore, ∇ 2 h S (w) ΛI from Assumption 1. Consequently, for every b ∈ R d we have Minimizing the right hand side of the above in b gives Re-arranging the above and switching back to f , we have Recalling the notation δ f (x) def = f (x) − f (w * ) and taking expectation with respect to S gives Now we apply the above to obtain where we used a + b 2 2 ≤ 2 a 2 2 + 2 b 2 2 in the first inequality and µ = E [∇f S (w k ) − ∇f S (w * )]. Finally, E g 2 2 (18) ≤ 4Λ (δ f (x) + δ f (w)) − 2 ∇f (w) 2 2 ≤ 4Λδ f (x) + 4(Λ − λ)δ f (w).
In the last inequality we used the fact that strongly convex functions satisfy ∇f (x) 2 2 ≥ 2λ δ f (x) for all x ∈ R d .
8 Proof of Lemma 1 To simplify notation, we define G Proposition 2.2 in [11] proves that so long as G and H (and hence B) are positive definite and D has full rank, then H + is positive definite and non-singular, and consequently, B + is well defined and positive definite. Using the Sherman-Morrison-Woodbury identity the update formula for B + , as shown in the Appendix in [11], is given by