Iterative Methods for Sylvester-like Variational Inequality Problems
Ning Zheng
Abstract
We consider the solution of Sylvester-like variational inequality problem (SVIP), or linear matrix equation complementarity problem (LMECP),
\(\seteqnumber{0}{}{0}\)\begin{align} \label {eqn:lmecp} X\geq \bzero ,\quad W=AX+XB-F\geq \bzero \quad \mbox {and}\quad \langle X, W\rangle = \bzero \end{align} where \(A\in {\bf R}^{m\times m}\), \(B\in {\bf R}^{n\times n}\) and \(F\in {\bf R}^{m\times n}\) are large, sparse and structured discretization matrices from partial differential operators, and \(X\in {\bf R}^{m\times n}\) is an unknown matrix. Here, \(\langle X,W \rangle ={\rm Tr}(X^{\top }W)\) denotes the Frobenius inner product of two matrices, where \(X^{\top }\) denotes the transpose of the matrix \(X\). If \(B=A^{\top }\) and \(F\) is symmetric, we refer (1) as the Lyapunov-like linear matrix equation complementarity problem. LMECP (1) generally arises from the finite discretization of free boundary problems
\(\seteqnumber{0}{}{1}\)\begin{align} \label {eqn:lcppde} v(x)\geq g(x),\quad w(x)=\mathcal {L}v(x)-f(x)\geq 0 \quad \mbox {and}\quad (v(x)-g(x))w(x)=0, \end{align} where \(\mathcal {L}\) is a given partial differential operator, and \(x\in D\subseteq {\bf R}^n\) where \(D\) is a given domain with boundary \(\partial D\). The boundary condition of (2) is \(v(x)=g(x),~x\in \partial D\). Well known examples of free boundary problems which can be written in the form (2) include American option pricing, porous flow through dams, journal bearing lubrication, and elastic-plastic torsion, etc.
The vectorization of the LMECP (1) gives a mathematically equivalent linear complementarity problem (LCP),
\(\seteqnumber{0}{}{2}\)\begin{align} \label {eqn:lmecp-vec} \bx \geq \bzero ,\quad \mathcal {A}\bx -\bof \geq \bzero \quad \mbox {and}\quad \bx ^{\top }(\mathcal {A}\bx -\bof )=0, \end{align} where \(\mathcal {A}=I_n\otimes A+B^{\top }\otimes I_m\in {\bf R}^{mn\times mn}\), \(\bof ={\rm vec}(F)\in {\bf R}^{mn\times 1}\) and \(\bx ={\rm vec}(X)\in {\bf R}^{mn\times 1}\). Here, the symbol \(\otimes \) denotes the Kronecker product, and \({\rm vec}(\cdot )\) denotes the vectorization operator that converts a matrix into a vector by stacking the columns of the matrix on top of one another. There are few numerical methods specifically designed for solving LMECP (1). Numerical methods [6, 4, 2, 3, 7] for LCP (3) are generally not efficient for directly solving LMECP (1) due to the large storage and complexity.
When the matrix \(\mathcal {A}\) arises from the finite difference discretization of elliptic and parabolic partial differential equations, it has structure that contains discretization components from different spatial derivatives. Hence, the idea of alternating direction implicit (ADI) method is to split the finite difference operator into separate operators, where each operator corresponds to the discretization of one-dimensional spatial derivative term, so that the solution of discretized system can be transformed to the alternative solutions of discretized sub-systems, which have a simpler structure that requires fewer storage and computational costs. Let \(\mathcal {A}=\mathcal {H}+\mathcal {V}\) be the matrix splitting of \(\mathcal {A}\), where \(\mathcal {H}=I_n\otimes A\) and \(\mathcal {V}=B^{\top }\otimes I_m\) are respectively generated from discrete central difference approximations to the particular one-dimensional equation. The Peaceman-Rachford ADI for LCP (3) proposed by Lin and Cryer [5], and the matricization of ADI is described as follows.
Algorithm 1: Peaceman-Rachford ADI for LMECP (1)
-
1: Input an initial guess \(X^{(0)}\) and positive parameters \(r_k\)
-
2: for \(k=0, 1, 2, \cdots \) until convergence do
-
3: Compute \(X^{\left (k+\frac {1}{2}\right )}\) by solving the LCP subproblem
\(\seteqnumber{0}{}{3}\)\begin{align} \left \{ \begin{array}{ll} W^{\left (k+\frac {1}{2}\right )}=(A+r_k I_m)X^{\left (k+\frac {1}{2}\right )}+X^{(k)}(B-r_k I_n)-F \geq \bzero ,& \\ X^{\left (k+\frac {1}{2}\right )}\geq \bzero ,\quad \left \langle X^{\left (k+\frac {1}{2}\right )}, W^{\left (k+\frac {1}{2}\right )}\right \rangle =0 & \end {array} \right . \label {eqn:ADI-LMECP-1} \end{align}
-
4: Compute \(X^{(k+1)}\) by solving the LCP subproblem
\(\seteqnumber{0}{}{4}\)\begin{align} \left \{ \begin{array}{ll} W^{(k+1)}=X^{(k+1)}(B+r_k I_n)+(A-r_k I_m)X^{\left (k+\frac {1}{2}\right )}-F \geq \bzero ,& \\ X^{(k+1)}\geq \bzero ,\quad \left \langle W^{(k+1)}, X^{(k+1)}\right \rangle =0. & \end {array} \right . \label {eqn:ADI-LMECP-2} \end{align}
-
5: end for
First, we propose a projection method for solving LMECP (1) by transforming the matrix equation into LCP (3) with a vector form by means of the Kronecker product. We can reformulate LCP (3) to an equivalent fixed-point equation
\(\seteqnumber{0}{}{5}\)\begin{align*} \bx ={\bf Proj}(\bx -\alpha [\mathcal {A}\bx -\bof ]), \end{align*} where \(\alpha >0\) is a scalar and \({\bf Proj}(\cdot )=\max (\cdot ,0)\) denotes the orthogonal projection of vector or matrix onto nonnegative cone. The matricization form gives
\(\seteqnumber{0}{}{5}\)\begin{align*} X={\bf Proj}(X-\alpha (AX+XB-F)). \end{align*}
The gradient projection method for LMECP (1) is listed in Algorithm 2.
Algorithm 2: Projection method for LMECP (1)
-
1: Input an initial guess \(X^{(0)}\) and positive parameter \(\alpha \)
-
2: for \(k=0, 1, 2, \cdots \) until convergence do
-
3: Compute the residual \(R^{(k)}=F-AX^{(k)}-X^{(k)}B\)
-
4: Compute
\(\seteqnumber{0}{}{5}\)\begin{align} \label {eqn:projection-iteration} X^{(k+1)}={\bf Proj}(X^{(k)}+\alpha R^{(k)}). \end{align}
-
5: end for
Next, we discuss the convergence of Peaceman-Rachford ADI algorithm Algorithm 1 for the non-Hermitian case. Unlike the symmetric case [1, 5], convergence properties for nonsymmetric situations cannot be established relying on the descent function of the quadratic form. Rather, as with most iterative methods for solving systems of equations, the recursive relation between two successive iterations will be utilized here. We first equivalently reformulate the LCP (1) as an implicit fixed-point equation by variable transformation, and thus the ADI Algorithm 1 can be correspondingly reformulated. Then the recursive error relations are constructed based on the fixed-point equations. In addition, we consider the case when \(\mathcal {H}\) and \(\mathcal {V}\) are \(H_+\)-matrices. We study the convergence analysis of ADI algorithm for LMECP when \(\mathcal {H}\mathcal {V}=\mathcal {V}\mathcal {H}\) does not necessary hold.
Denote
\(\seteqnumber{0}{}{6}\)\begin{align*} L_{\alpha }(\mathcal {H}) &=|(\alpha I+\mathcal {H}+r_kI)^{-1}(\alpha I-\mathcal {H}-r_kI)|, \\ L_{\beta }(\mathcal {V}) &=|(\beta I+\mathcal {V}+r_kI)^{-1}(\beta I-\mathcal {V}-r_kI)|, \\ K_{\alpha }(\mathcal {H},\mathcal {V}) &=2|(\alpha I+\mathcal {H}+r_kI)^{-1}(\mathcal {V}-r_k I)| \\ K_{\beta }(\mathcal {V},\mathcal {H}) &=2|(\beta I+\mathcal {V}+r_kI)^{-1}(\mathcal {H}-r_k I)| \end{align*} We have the following convergence theorem.
-
Theorem 1 Peaceman-Rachford ADI algorithm converges to the unique solution for any initial vector if \(\rho (L_{\alpha }(\mathcal {H}))<1\), \(\rho (L_{\beta }(\mathcal {V}))<1\) and \(\rho (G)<1\), where \(\rho (\cdot )\) denotes for the spectral radius of the matrix and
\(\seteqnumber{0}{}{6}\)\begin{align*} G=[I-L_{\beta }(\mathcal {V})]^{-1}K_{\beta }(\mathcal {V},\mathcal {H})[I-L_{\alpha }(\mathcal {H})]^{-1}K_{\alpha }(\mathcal {H},\mathcal {V}). \end{align*}
Consider the case when both \(\mathcal {H}\) and \(\mathcal {V}\) are \(H_+\)-matrices. Let \(\mathcal {H}=D_{\mathcal {H}}+B_{\mathcal {H}}\) and \(\mathcal {V}=D_{\mathcal {V}}+B_{\mathcal {V}}\), where \(D_\mathcal {H}\) and \(B_\mathcal {H}\) are the diagonal and off-diagonal parts of \(\mathcal {H}\), respectively, and \(D_\mathcal {V}\) and \(B_\mathcal {V}\) are the diagonal and off-diagonal parts of \(V\), respectively.
-
Theorem 2 Peaceman-Rachford ADI algorithm converges to the unique solution for any initial vector, provided that \(\mathcal {H}\) and \(\mathcal {V}\) are \(H_+\)-matrices and
\(\seteqnumber{0}{}{6}\)\begin{align*} (\alpha -r_k)I-D_\mathcal {H}\geq \bzero & \quad \mbox {and}\quad D_\mathcal {V}-r_kI\leq \bzero , \\ (\beta -r_k)I-D_\mathcal {V}\geq \bzero & \quad \mbox {and}\quad D_\mathcal {H}-r_kI\leq \bzero , \end{align*}
In the following analysis, we assume that \(\mathcal {H}\) and \(\mathcal {V}\) are commute, that is to say, \(\mathcal {H}\mathcal {V}=\mathcal {V}\mathcal {H}\). Remark that for \(\mathcal {H}\) and \(\mathcal {V}\) arising from the finite difference discretization of a separable second-order elliptic operator in a rectangular region, it can be shown that \(\mathcal {H}\mathcal {V}=\mathcal {V}\mathcal {H}\) holds.
Instead of taking the absolute value, we give another general convergence result based on the matrix norm as follows. Denote
\(\seteqnumber{0}{}{6}\)\begin{align*} \delta _{\alpha }(\mathcal {H}) & = \|(\alpha I+\mathcal {H}+r_kI)^{-1}(\alpha I-\mathcal {H}-r_kI)\|, \\ \delta _{\beta }(\mathcal {V}) & = \|(\beta I+\mathcal {V}+r_kI)^{-1}(\beta I-\mathcal {V}-r_kI)\|, \\ \tau _{\alpha }(\mathcal {H}) & =2\|(\alpha I+\mathcal {H}+r_kI)^{-1}(\mathcal {V}-r_k I)\|, \\ \tau _{\beta }(\mathcal {V}) & =2\|(\beta I+\mathcal {V}+r_kI)^{-1}(\mathcal {H}-r_k I)\|, \end{align*} where \(\|\cdot \|\) denotes for matrix norm.
-
Theorem 3 Peaceman-Rachford ADI algorithm converges to the unique solution for any initial vector if
\(\seteqnumber{0}{}{6}\)\begin{align} \label {eqn:convcondition} \delta _{\alpha }(\mathcal {H})<1,\quad \delta _{\beta }(\mathcal {V})<1 \quad \mbox {and}\quad \frac {\tau _{\alpha }(\mathcal {H})\tau _{\beta }(\mathcal {V})}{[1-\delta _{\alpha }(\mathcal {H})][1-\delta _{\beta }(\mathcal {V})]}<1. \end{align}
Consider the case when both \(\mathcal {H}\) and \(\mathcal {V}\) are Hermitian positive definite matrices, and thus \(\mathcal {A}=\mathcal {H}+\mathcal {V}\) is Hermitian positive definite.
-
Theorem 4 Suppose that \(\mathcal {H}\) and \(\mathcal {V}\) are Hermitian positive definite matrices, and \(\mathcal {H}\) and \(\mathcal {V}\) are commute. If
\(\seteqnumber{0}{}{7}\)\begin{equation*} r_k\geq \frac {1}{2}\max (\lambda _1+\lambda _n,\sigma _1+\sigma _n), \end{equation*}
then Peaceman-Rachford ADI algorithm for LCP converges to the unique solution for any initial vector.
Finally, we present numerical experiments to show the convergence of proposed methods. We consider the free boundary problem arises from fractional Black-Scholes American option pricing. Assume that the asset prices \(S_1\) and \(S_2\) satisfy independent Lévy stochastic processes
\(\seteqnumber{0}{}{7}\)\begin{align*} \mathcal {L}u=-\frac {\partial u}{\partial t}+a_1 \frac {\partial u}{\partial x}+a_2 \frac {\partial u}{\partial y}-b_1\left [ \prescript {}{-\infty }{D_x^\alpha u}\right ]-b_2\left [\prescript {}{-\infty }{D_y^\beta u}\right ]+r u \end{align*} where \(x=\ln S_1\) and \(y=\ln S_2\). Here, \(\prescript {}{-\infty }{D_x^\alpha u}\) and \(\prescript {}{-\infty }{D_y^\beta u}\) represents Caputo derivatives of \(u\) on \(x\) and \(y\), and \(\alpha , \beta \in (1,2)\).
\[ \begin {array}{ll} a_1=-r-\frac {1}{2} \sigma _1^\alpha \sec \left (\frac {\alpha \pi }{2}\right ), & b_1=-\frac {1}{2} \sigma _1^\alpha \sec \left (\frac {\alpha \pi }{2}\right ) \\ a_2=-r-\frac {1}{2} \sigma _2^\beta \sec \left (\frac {\beta \pi }{2}\right ), & b_2=-\frac {1}{2} \sigma _2^\beta \sec \left (\frac {\beta \pi }{2}\right ) \end {array} \]
\(\sigma _1, \sigma _2\) are volatilities of asset prices. By finite difference discretization of the model, we apply the projection method and ADI algorithm to solve the resulting LMECP, and the numerical results further confirm our convergence analysis.
References
-
[1] B. H. Ahn, Solution of nonsymmetric linear complementarity problems by iterative methods, Journal of Optimization Theory and Applications, 33 (1981), pp. 175–185.
-
[2] Z.-Z. Bai, Modulus-based matrix splitting iteration methods for linear complementarity problems, Numerical Linear Algebra with Applications, 6 (2010), pp. 917–933.
-
[3] A. Hadjidimos and M. Tzoumas, Nonstationary extrapolated modulus algorithms for the solution of the linear complementarity problem, Linear Algebra and Its Applications, 12 (2009), pp. 197–210.
-
[4] N. W. Kapple and L. T. Watson, Iterative algorithms for the linear complementarity problem, International Journal of Computer Mathematics, 19 (1986), pp. 273–297.
-
[5] Y. Lin and C. W. Cryer, An alternating direction implicit algorithm for the solution of linear complementarity problems arising from free boundary problems, Applied Mathematics and Optimization, 13 (1985), pp. 1–17.
-
[6] K. G. Murty and F.-T. Yu, Linear complementarity, linear and nonlinear programming, volume 3, Citeseer, 1988.
-
[7] N. Zheng, K. Hayami and J.-F. Yin, Modulus-type inner outer iteration methods for nonnegative constrained least squares problems, SIAM Journal on Matrix Analysis and Applications, 37 (2016), pp. 1250–1278.