A framework of constraint preserving update schemes for optimization on Stiefel manifold

Abstract

This paper considers optimization problems on the Stiefel manifold $X^{\mathsf{T}}X=I_p$, where $X\in \mathbb {R}^{n \times p}$ is the variable and $I_p$ is the $p$-by-$p$ identity matrix. A framework of constraint preserving update schemes is proposed by decomposing each feasible point into the range space of $X$ and the null space of $X^{\mathsf{T}}$. While this general framework can unify many existing schemes, a new update scheme with low complexity cost is also discovered. Then we study a feasible Barzilai–Borwein-like method under the new update scheme. The global convergence of the method is established with an adaptive nonmonotone line search. The numerical tests on the nearest low-rank correlation matrix problem, the Kohn–Sham total energy minimization and a specific problem from statistics demonstrate the efficiency of the new method. In particular, the new method performs remarkably well for the nearest low-rank correlation matrix problem in terms of speed and solution quality and is considerably competitive with the widely used SCF iteration for the Kohn–Sham total energy minimization.

Appendices

Appendix 1: Details of approach I and II in Sect. (3)

1.1 Details of approach I

Using the condition that $Z(\tau ) \equiv I_p$ and $\widetilde{R}'(0) = - X^{\mathsf{T}} E$ and denoting ${\mathcal {A}}(\tau ) = \left[ \begin{array}{c} \widetilde{R}(\tau ) \\ W \widetilde{N}(\tau ) \end{array}\right]$, it follows from (3.4) and (3.6) that

$$\begin{aligned} {\mathcal {A}}'(0) = \left[ \begin{array}{cc} -X^{\mathsf{T}}E &{}~B \\ W &{} ~F \end{array}\right] \left[ \begin{array}{c} I_p \\ 0\end{array}\right] \ \mathrm {and} \ {\mathcal {A}}(0) = \left[ \begin{array}{c} I_p \\ 0\end{array}\right] , \end{aligned}$$

where $B \in {\mathbb {R}}^{n\times p}$, $F \in {\mathbb {R}}^{p \times p}$. Solving the above ordinary differential equations, we get that

$$\begin{aligned} {\mathcal {A}}(\tau ) = \exp \left( \tau \left[ \begin{array}{cc} -X^{\mathsf{T}}E &{}~B \\ W &{} ~F \end{array}\right] \right) \left[ \begin{array}{c} I_p \\ 0\end{array}\right] . \end{aligned}$$

(8.1)

Since $Z(\tau ) \equiv I_p$, we know from (3.3) and the definition of ${\mathcal {A}}(\tau )$ that

$$\begin{aligned} Y(\tau ) = \big [\begin{array}{l} X, \ I_p\end{array} \big ] {\mathcal {A}}(\tau ), \end{aligned}$$

and ${\mathcal {A}}(\tau )^{\mathsf{T}} {\mathcal {A}}(\tau ) = I_p$. It follows that the matrix inside the brackets of (8.1) will be skew-symmetric. This means that $B = -W^{\mathsf{T}}$, $F + F^{\mathsf{T}} = 0$ and $X^{\mathsf{T}}E+E^{\mathsf{T}}X=0$ (i.e., $E \in {\mathcal {T}}_X$). Let $W = QR$ be the unique QR factorization of $W$ and assume that $F$ takes the form $F = Q{\widetilde{F}}Q^{\mathsf{T}}$, where ${\widetilde{F}}\in {\mathbb {R}}^{p\times p}$ is any skew-symmetric matrix. Then we can write

$$\begin{aligned} Y(\tau ) =&\big [\begin{array}{l} X, \ I_p\end{array} \big ] \exp \left( \tau \left[ \begin{array}{cc} I_p &{}~0 \\ 0 &{} ~Q \end{array}\right] \left[ \begin{array}{cc} -X^{\mathsf{T}}E &{}~-R^{\mathsf{T}} \\ R &{} ~{\widetilde{F}}\end{array}\right] \left[ \begin{array}{cc} I_p &{}~0 \\ 0 &{} ~Q^{\mathsf{T}} \end{array}\right] \right) \left[ \begin{array}{c} I_p \\ 0\end{array}\right] \nonumber \\ =&\big [\begin{array}{l} X, \ Q \end{array} \big ] \exp \left( \tau \left[ \begin{array}{cc} -X^{\mathsf{T}}E &{}~-R^{\mathsf{T}} \\ R &{} ~{\widetilde{F}}\end{array}\right] \right) \left[ \begin{array}{c} I_p \\ 0\end{array}\right] , \end{aligned}$$

(8.2)

where the second equality is due that $Q$ has orthonormal columns. The update scheme (8.2) can be regarded as a generalized geodesic update scheme, since letting ${\widetilde{F}}= 0$, (8.2) reduces to the geodesic update scheme (2.1).

1.2 Details of approach II

Based on the choices that $\widetilde{R}(\tau ) = I_p + \tau \widetilde{R}'(0)$ and $\widetilde{N}(\tau ) = \tau I_p$, we can get $Z(\tau )$ from (3.3) by the polar decomposition or the Cholesky factorization.

If by the polar decomposition, $Z(\tau )$ and $Z'(\tau )$ are always symmetric. In this case, it is easy to show that $Y(\tau ; X) = \mathcal {P}_{{\mathcal {S}}_{n,p}}{{\big (X \widetilde{R}(\tau ) + W \widetilde{N}(\tau ) \big )}}$, which with (3.6), (3.7) and (3.8) further implies that

$$\begin{aligned} Y(\tau ; X) = \mathcal {P}_{{\mathcal {S}}_{n,p}}\big (X - \tau E - \tau X Z'(0) \big ), \end{aligned}$$

(8.3)

where $Z'(0)$ is any $p$-by-$p$ symmetric matrix. If $Z'(0) = 0$, (8.3) reduces to the ordinary polar decomposition or Manton’s projection update scheme. If $Z'(0) = \text{ sym }(X^{\mathsf{T}}G)$ and $E = G - X \text{ sym }(X^{\mathsf{T}}G),$ it becomes the ordinary gradient projection update scheme.

If by the Cholesky factorization, $Z(\tau )$ and $Z'(\tau )$ are always upper triangular. Similarly, we can derive

$$\begin{aligned} Y(\tau ; X) = \text{ qr }\big (X - \tau E - \tau X Z'(0)\big ), \end{aligned}$$

(8.4)

where $Z'(0)$ is any $p$-by-$p$ upper triangular matrix. When $Z'(0) = 0$, (8.4) reduces to the ordinary QR factorization update scheme.

Appendix 2: Proof of lemma 3.1

The fact that $X^{\mathsf{T}}D_{\rho }$ is skew-symmetric implies that $J$ is invertible, thus $Y(\tau )$ is well-defined. (i) follows from the construction of the update scheme (3.11).

Meanwhile, we know that $Y'(0) = -D_{\rho }$, which with the chain rule shows that

$$\begin{aligned} {\mathcal {F}}'_{\tau }(Y(0)) = - \langle G, D_{\rho }\rangle = -\langle \nabla {\mathcal {F}}, P_X D_{\rho }\rangle , \end{aligned}$$

where the second equality is due to $D_{\rho }\in {\mathcal {T}}_X$ and the definition (2.8) of $\nabla {\mathcal {F}}$. Recall that $P_X = I_n - \frac{1}{2} XX^{\mathsf{T}}$. Substituting (3.12) into the above equation yields

$$\begin{aligned} {\mathcal {F}}'_{\tau }(Y(0))= -\langle \nabla \mathcal {F}, (I_n + (\rho -1 )XX^{\mathsf{T}})\nabla \mathcal {F}\rangle \le - \min \{\rho ,1\} \Vert \nabla \mathcal {F}\Vert _{\mathsf{F}}^2. \end{aligned}$$

So (ii) is true.

We prove (iii) by contradiction. Multiplying $X^{\mathsf{T}}$ from both sides the expression of $Y(\tau )$ in (3.11) and using $X^{\mathsf{T}}W = 0$, we get that $2J^{-1} - I_p = X^{\mathsf{T}}Y$. Assume that there exists a $p$-by-$p$ orthogonal matrix $Q_p$ such that $Y = XQ_p$. Then we have $2J^{-1} - I_p = Q_p$; i.e., $2I_p - J = Q_pJ$. It follows from $(2I_p - J)^{\mathsf{T}} (2I_p-J) = (Q_pJ)^{\mathsf{T}}Q_pJ$ and $Q_p^\mathsf{T}Q_p = I_p$ that

$$\begin{aligned} J^{\mathsf{T}} + J = 4I_p, \end{aligned}$$

which is a contradiction to the definition of $J$. The contradiction shows the truth of (iii).

For (iv), it is obvious that $\Vert J\Vert _2 \le 1 + \frac{\tau ^2}{4}\Vert W\Vert _{\mathsf{F}}^2 + \frac{\tau }{2}\Vert X^{\mathsf{T}}D_{\rho }\Vert _{\mathsf{F}}$, where $\Vert \cdot \Vert _2$ is the spectral norm. It follows from $-D_{\rho }= W - XX^{\mathsf{T}} D_{\rho }$ and $X^{\mathsf{T}} W=0$ that $\Vert D_{\rho }\Vert _{\mathsf{F}}^2 = \Vert W\Vert _{\mathsf{F}}^2 + \Vert X^{\mathsf{T}}D_{\rho }\Vert _{\mathsf{F}}^2.$ Notice that $\upsilon = \tau \Vert D_{\rho }\Vert _{\mathsf{F}}$. Then we have

$$\begin{aligned} \Vert J\Vert _2&\le 1 + \frac{\upsilon ^2}{4} \cdot \frac{\Vert W\Vert _{\mathsf{F}}^2}{\Vert D_{\rho }\Vert _{\mathsf{F}}^2 } + \frac{\upsilon }{2} \cdot \frac{\Vert X^{\mathsf{T}} D_{\rho }\Vert _{\mathsf{F}}^{}}{\Vert D_{\rho }\Vert _{\mathsf{F}}^{}}\\&\le 1 + \max _{0 \le t \le 1}\left( \frac{\upsilon ^2}{4} (1 - t^2) + \frac{\upsilon }{2} t\right) \le (5 + \upsilon ^2)/4. \end{aligned}$$

On the other hand, the relation $2J^{-1} = X^{\mathsf{T}} Y + I_p$ indicates that $\Vert J^{-1}\Vert _2 \le 1$. Thus (iv) is true.

In the case that $p=1$ or $n$, it is not difficult to simplify the corresponding update schemes and we omit the details here. Thus, we complete the proof.

Appendix 3: Proof of lemma 3.2

The relation $\langle G, D^{(i)} \rangle = e_i^{\mathsf{T}}\left( G^{\mathsf{T}} \nabla {\mathcal {F}}\right) e_i$ can be easily verified from the definition of $D^{(i)}$. With the definition of $D^{(q)}$ and $X^{\mathsf{T}} X = I_p$, we have that

$$\begin{aligned} {{X^{\mathsf{T}} D^{(q)} = X^{\mathsf{T}} G_{(q)}^{} e_q^{\mathsf{T}} - e_q G_{(q)}^{\mathsf{T}}X = \frac{2}{\tau }\left( be_q^{\mathsf{T}} - e_q b^{\mathsf{T}}\right) \in {\mathcal {T}}_X}}. \end{aligned}$$

Thus $Y(\tau )$ is well-defined. Moreover, we have that

$$\begin{aligned} W = -(I_n - XX^{\mathsf{T}}) D^{(q)} = -(I_n - XX^{\mathsf{T}}) G_{(q)}^{}e_q^{\mathsf{T}} \end{aligned}$$

Hence, the matrix $J = I_p + \frac{\tau ^2}{4} W^{\mathsf{T}}W + \frac{\tau }{2}X^{\mathsf{T}} D^{(q)}$ can be expressed as

$$\begin{aligned} J = I_p + \xi e_q^{}e_q^{\mathsf{T}} + b e_q^{\mathsf{T}} - e_q b^{\mathsf{T}} = {{I_p + \big [\begin{array}{cc} e_q,&\ b \end{array} \big ] \left( \left[ \begin{array}{cc} \xi &{}\ -1 \\ 1 &{} \ 0\end{array} \right] \left[ \begin{array}{c}e_q^{\mathsf{T}}\\ b^{\mathsf{T}}\end{array}\right] \right) }}, \end{aligned}$$

(10.1)

where $\xi = \frac{\tau ^2}{4}G_{(q)}^{\mathsf{T}}(I_n - XX^{\mathsf{T}}) G_{(q)}$. By the formulations of $b$ and $\xi$, we easily see that

$$\begin{aligned} {{e_q^{\mathsf{T}} b = 0, \ \xi =\alpha \, - \, b^{\mathsf{T}}b.}} \end{aligned}$$

Applying the SMW formula to (10.1) and using the above relations, we can obtain (3.15).

Moreover, notice that $\nabla {\mathcal {F}}= \sum _{i=1}^{p} D^{(i)}$. Thus we have

$$\begin{aligned} p\, {\mathcal {F}}'_{\tau }(Y(0)) = -p\, \langle G, D^{(q)} \rangle \le - \langle G, \nabla {\mathcal {F}}\rangle \le -\frac{1}{2}\Vert \nabla {\mathcal {F}}\Vert _{\mathsf{F}}^2, \end{aligned}$$

where the first inequality follows from the choice of $q$ in (3.14) and the second one is due to (ii) of Lemma 3.1. The proof is completed.

Appendix 4: Proof of proposition 3.1

Before going into the details of the proof, we recall the following fact on the inverse of a $2 \times 2$ block matrix.

Fact 11.1

Assume that $T = \left[ \begin{array}{cc} T_{11} &{} ~T_{12}\\ T_{21} &{}~T_{22} \end{array} \right]$, where $T_{22}^{}$ and $T_{11}^{} - T_{12}^{} T_{22}^{-1}T_{21}^{}$ are invertible. Then $T$ is invertible and

$$\begin{aligned} T^{-1} =\left[ \begin{array}{ll} ( T_{11}^{} - T_{12}^{} T_{22}^{-1}T_{21}^{})^{-1} &{} - ( T_{11}^{} - T_{12}^{} T_{22}^{-1}T_{21}^{})^{-1} T_{12}^{}T_{22}^{-1} \\ -T_{22}^{-1}T_{21}^{} ( T_{11}^{} - T_{12}^{} T_{22}^{-1}T_{21}^{})^{-1} &{}~ ~ T_{22}^{-1}T_{21}^{} ( T_{11}^{} - T_{12}^{} T_{22}^{-1}T_{21}^{})^{-1} T_{12}^{}T_{22}^{-1} + T_{22}^{-1}\end{array} \right] . \end{aligned}$$

First, we show that the update scheme (2.3) is well-defined, provided that $I_p+\frac{\tau }{4}X^{\mathsf{T}}D$ is invertible. Consider the update scheme (2.3) with $U=[P_X D, \,X]$ and $V=[X,\, -P_XD]$. It follows from $P_X = I_n - \frac{1}{2}XX^{\mathsf{T}}$ that $X^{\mathsf{T}}P_X D = \frac{1}{2} X^{\mathsf{T}}D$. Combining this with $X^{\mathsf{T}}X = I_p$ and $X^{\mathsf{T}}D$ being skew-symmetric, we can rewrite

$$\begin{aligned} I_{2p}+\frac{\tau }{2}V^{\mathsf{T}}U=\left[ \begin{array}{cc}I_p+\frac{\tau }{4}X^{\mathsf{T}}D &{} \frac{\tau }{2}I_p\\ -\frac{\tau }{2}D^{\mathsf{T}}P_X^{\mathsf{T}}P_X^{}D &{}~ I_p +\frac{\tau }{4}X^{\mathsf{T}}D\end{array}\right] . \end{aligned}$$

(11.1)

Moreover, we derive that

$$\begin{aligned} {{W^{\mathsf{T}}W = D^{\mathsf{T}}(I_n - XX^{\mathsf{T}})D = D^{\mathsf{T}}P_X^{\mathsf{T}} P_X^{} D + \frac{1}{4} \left( X^{\mathsf{T}}D\right) ^2,}} \end{aligned}$$

which with (3.11) implies that

$$\begin{aligned} J = \Big (I_p + \frac{\tau }{4}X^{\mathsf{T}} D \Big )^2 + \frac{\tau ^2}{4}D^{\mathsf{T}}P_X^{\mathsf{T}}P_X^{} D. \end{aligned}$$

Plugging this into (11.1) yields

$$\begin{aligned} I_{2p}+\frac{\tau }{2}V^{\mathsf{T}}U:= \left[ \begin{array}{cc}T_{11} &{} ~T_{12} \\ T_{21} &{}~ T_{22} \end{array}\right] =\left[ \begin{array}{cc}I_p+\frac{\tau }{4}X^{\mathsf{T}}D &{} \frac{\tau }{2}I_p\\ \frac{2}{\tau } \left( \left( I_p +\frac{\tau }{4}X^{\mathsf{T}}D\right) ^2 - J\right) &{}~ I_p +\frac{\tau }{4}X^{\mathsf{T}}D\end{array}\right] , \end{aligned}$$

(11.2)

where $T_{ij} \in {\mathbb {R}}^{p\times p}$ $(i,j \in \{1,2\})$. By simple calculations, we know that

$$\begin{aligned} {{T_{11}^{} - T_{12}^{} T_{22}^{-1} T_{21}^{} = \left( I_p + \frac{\tau }{4}X^{\mathsf{T}}D\right) ^{-1} J}}. \end{aligned}$$

Thus it follows from Fact 11.1 that $I_{2p}+\frac{\tau }{2}V^{\mathsf{T}}U$ is invertible and the update scheme (2.3) is well-defined.

We now prove the equivalence. With Lemma 4.1 and (11.2), some tedious manipulations yield that $(I_{2p}+\frac{\tau }{2}V^{\mathsf{T}}U)^{-1}=\left[ \begin{array}{cc}M_{11} &{} ~M_{12}\\ M_{21} &{} ~M_{22}\end{array}\right]$, where

$$\begin{aligned}\begin{array}{ll} M_{11} = J^{-1}(I_p + \frac{\tau }{4}X^{\mathsf{T}} D), &{} ~ M_{12} = -\frac{\tau }{2}J^{-1}, \\ M_{21} = \frac{2}{\tau }\big (I_p- {{M_{22}}}(I_p+\frac{\tau }{4}X^{\mathsf{T}}D)\big ), &{}~ M_{22}=(I_p+\frac{\tau }{4}X^{\mathsf{T}}D)J^{-1}. \end{array} \end{aligned}$$

Direct calculations show that

$$\begin{aligned} M_{11}+ \frac{1}{2}M_{12}X^{\mathsf{T}}D = J^{-1}, \quad \tau \Big (M_{21} + \frac{1}{2} M_{22}X^{\mathsf{T}}D\Big ) = 2I_p - 2\Big (I_p+ \frac{\tau }{4}X^{\mathsf{T}}D\Big )J^{-1}. \end{aligned}$$

(11.3)

Finally, we can obtain

$$\begin{aligned} Y_{\mathrm {wy}}(\tau ) ={}&X - \tau \left[ \begin{array}{l}P_XD,\ X\end{array}\right] \Bigg (\left[ \begin{array}{cc}M_{11} &{} ~M_{12}\\ M_{21} &{} ~M_{22}\end{array}\right] \left[ \begin{array}{c} I_p\\ \frac{1}{2}X^{\mathsf{T}}D \end{array}\right] \Bigg ) \nonumber \\ ={}&X-\tau \Big (\frac{1}{2} XX^{\mathsf{T}}D -W\Big )\Big (M_{11}+ \frac{1}{2}M_{12}X^{\mathsf{T}}D\Big ) -\tau X\Big (M_{21}+\frac{1}{2}M_{22}X^{\mathsf{T}}D\Big )\nonumber \\ ={}&(2X + \tau W)J^{-1}- X, \end{aligned}$$

where the first equality uses (2.3) and $X^{\mathsf{T}}X = I_p$, the second and third equalities are due to $P_XD = \frac{1}{2} XX^{\mathsf{T}}D -W$ and (11.3), respectively. The proof is completed.

Appendix 5: Some details of Table 1

We first review the computational costs of some basic matrix operations. Given $A \in {\mathbb {R}}^{n_1 \times n_2}$ and $B \in {\mathbb {R}}^{n_2 \times n_2}$, computing $AB$ needs $2n_1^{}n_2^2$ flops while computing $A^{\mathsf{T}} A$ only needs $n_1^{}n_2^2$ flops. Computing $\text{ qr }(A)$ by the modified Gram-Schmidt algorithm needs $2n_1^{}n_2^2$ flops (here $n_1\ge n_2$). If $B$ is symmetric, computing the eigenvalue decomposition $B = P\Sigma P^{\mathsf{T}}$ with orthogonal $P\in {\mathbb {R}}^{n_2 \times n_2}$ and diagonal $\Sigma \in {\mathbb {R}}^{n_2 \times n_2}$ by the symmetric QR algorithm needs $9n_2^3$ flops. If $B$ is skew-symmetric, computing the exponential of $B$ needs $10n_2^3$ flops. If $B$ is symmetric positive definite, computing $B^{1/2}$ needs $10n_2^3$ flops. In addition, if $B$ is nonsingular, solving the matrix equation $B T = A^{\mathsf{T}}$ by the Gauss elimination with partial pivoting needs $2n_1^{}n_2^2 +2n_2^3/3$ flops. It follows that computing $B^{-1}$ needs $8n_2^3/3$ flops. We refer interested readers to [23] for more details.

To verify the computational costs for the corresponding update scheme in Table 1, we take the new update scheme (3.11) and the Wen-Yin update scheme (2.3) as examples. Notice that comparing with $O(np^2)$ or $O(n^3)$, the $O(p^2)$ term is omitted for $1\le p \le n$. For simplicity, we only consider the case that $1< p <n$. In the case that $p=1$ or $n$, the cost for the update schemes (3.11) and (2.3) can be easily obtained by a similar analysis. We omit its details here.

To analyze the computational cost for (3.11), since $E = D_{\rho }$ and hence $W = -G + XX^{\mathsf{T}}G$, we can rewrite the feasible curve as

$$\begin{aligned} Y(\tau ; X) = \left( X\Big (\frac{2}{\tau }I_p + X^{\mathsf{T}}G\Big ) - G\right) {\left( \frac{J(\tau )}{\tau }\right) ^{-1}} - X, \end{aligned}$$

and

$$\begin{aligned} J(\tau ) = I_p + \frac{\tau ^2}{4} (G^{\mathsf{T}}G - G^{\mathsf{T}}XX^{\mathsf{T}}G) + {{\rho \tau }}(X^{\mathsf{T}}G - G^{\mathsf{T}}X). \end{aligned}$$

Forming $J(\tau )$ needs $3np^2 + p^3$ flops involving computing $X^{\mathsf{T}}G$, $G^{\mathsf{T}}G$ and $(G^{\mathsf{T}}X )(X^{\mathsf{T}}G )$. The final assembly for $Y(\tau )$ consists of involving solving one matrix equation and performing one matrix multiplication and two matrix subtractions, and hence needs $4np^2 + 2np + 2p^3/3$ flops. Therefore, the total cost of forming $Y(\tau ; X)$ in (3.11) for any $\tau$ is $7np^2 + 2np + 5p^3/3$. While doing backtracking line searches, we need to update $Y(\tau ; X)$ for a different $\tau$. Denote the first trial and the new trial stepsizes by $\tau _{{{\mathrm {first}}}}$ and $\tau _{{{\mathrm {new}}}}$, respectively. It is easy to see that

$$\begin{aligned} Y(\tau _{{{\mathrm {new}}}}; X) = \left( X\Big (\frac{2}{\tau _{{{\mathrm {first}}}}}I_p + X^{\mathsf{T}}G\Big ) - G + \Big (\frac{2}{\tau _{{{\mathrm {new}}}}} -\frac{2}{\tau _{{{\mathrm {first}}}}}\Big )X\right) {{\left( \frac{J(\tau _{{{\mathrm {new}}}})}{\tau _{\mathrm {new}}}\right) ^{-1}}} - X. \end{aligned}$$

As $X\big (\frac{2}{\tau _{{{\mathrm {first}}}}}I_p + X^{\mathsf{T}}G\big ) - G$ is already computed in $Y(\tau _{{{\mathrm {first}}}}; X)$, we store this matrix and hence only need $2np^2 + 3np + 2p^3/3$ to compute $Y(\tau _{{{\mathrm {new}}}}; X).$

To analyze the computational cost for the Wen-Yin update scheme (2.3), we see that it takes $3np^2$ and $2np^2 + np$ flops to form $I_{2p} + \frac{\tau }{2} V^{\mathsf{T}} U$ and $P_X D_{\rho }$, respectively. The final assembly for $Y(\tau )$ needs $4np^2 + np + 40p^3/3$ flops. Hence, the total cost for the Wen-Yin update scheme is $9np^2 + 2np + 40p^3/3$. When $\rho = 1/2,$ we see from [53] that $U = [G, X]$, $V = [X, -G]$ in (2.3) which implies that there is no need to compute $P_XD_{1/2}$ any more, and the total cost for forming the Wen-Yin update scheme can be reduced to $7np^2 + np + 40p^3/3$. The cost for updating $Y(\tau )$ with a new $\tau$ is $4np^2 + np + 40p^3/3$ since $V^{\mathsf{T}}U$ can be stored and the main work involves solving a matrix equation and the final assembly.