东华大学理学院，上海

收稿日期：2021年12月18日；录用日期：2022年1月20日；发布日期：2022年1月27日

摘要

线性矩阵不等式半定互补问题是一类新颖的多项式优化问题，是半定规划与互补问题的交叉研究内容。该问题可以转化为一系列带有线性矩阵不等式约束的多项式优化子问题，进而可以采用标量型松弛方法或矩阵型松弛方法进行求解。更进一步，当线性矩阵不等式半定互补问题的实数解个数有限时，利用本文给出的算法可以计算出该问题的全部实数解。最后，我们进行相关数值实验，分别使用标量型松弛方法和矩阵型松弛方法求解该问题，并将这两种方法的结果进行对比。

关键词

半定互补问题，线性矩阵不等式，标量型松弛方法，矩阵型松弛方法

Numerical Methods for the Semidefinite Complementarity Problem with Linear Matrix Inequalities

Shanshan Kong, Xin Zhao^*

College of Science, Donghua University, Shanghai

Received: Dec. 18^th, 2021; accepted: Jan. 20^th, 2022; published: Jan. 27^th, 2022

ABSTRACT

The semidefinite complementarity problem with linear matrix inequalities is a novel polynomial optimization problem, which is the cross-study content of semidefinite programming and complementarity problems. The problem can be transformed into a series of polynomial optimization subproblems with linear matrix inequality constraints. Then the problem can be solved by the scalar-type relaxation method or the matrix-type relaxation method. Furthermore, when the number of real solutions of the semidefinite complementarity problem with linear matrix inequalities is limited, all real solutions of the problem can be calculated by using the algorithm given in this paper. Finally, we conduct related numerical experiments. We use the scalar-type relaxation method and the matrix-type relaxation method to solve the problem and compare the results of these two methods.

Keywords:Semidefinite Complementarity Problem, Linear Matrix Inequalities, Scalar-Type Relaxation Method, Matrix-Type Relaxation Method

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

线性矩阵不等式(linear matrix inequalities, LMI)一般描述为：

$F\left(x\right)=F_0+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_i}{F}_i\succcurlyeq 0,$

其中 $x\in {ℝ}^n$ ， $x=\left(x_1,x_2,\cdots ,x_n\right)$ ， $F_i=F_i^{\text{T}}\in {ℝ}^{m\times m}$ ， $i=0,1,\cdots ,n$ 是对称矩阵，矩阵 $F\left(x\right)$ 是一个从n维实数空间到 $m\times m$ 维实对称矩阵空间的映射，矩阵中的每一个元素都是以 $x\in {ℝ}^n$ 为变量的多项式。矩阵 $F\left(x\right)\succcurlyeq 0$ 表示矩阵 $F\left(x\right)$ 是正定的，即对任意的 $u\in {ℝ}^m$ ，有 $u^{\text{T}}F\left(x\right)u\ge 0$ 。近年来，线性矩阵不等式得到了越来越多的关注和重视，被广泛地应用于解决控制中的问题，已经成为系统与控制领域研究中的一大热点研究内容 [1] [2]。

本文考虑带有线性矩阵不等式约束的半定互补问题：寻求解 $x\in {ℝ}^n$ 满足 $F_1\left(x\right)\succcurlyeq 0$ ， $F_2\left(x\right)\succcurlyeq 0$ ， $\langle F_1\left(x\right),F_2\left(x\right)\rangle =0$ ，其中 $F_i\left(x\right)\succcurlyeq 0,i=1,2$ 是线性矩阵不等式。线性矩阵不等式半定互补问题是一种新的半定互补问题，目前很少有关于这个问题的相关文献和理论结果。线性矩阵不等式半定互补问题可以看作标量型多项式互补问题的推广，这两个问题显著的区别是在线性矩阵不等式半定互补问题中涉及到矩阵型约束，而不仅仅是标量型约束。这里所说的“矩阵型”是指多项式矩阵不等式约束，“标量型”是指单个多项式不等式约束。

求解带有 $F\left(x\right)\succcurlyeq 0$ 的半定约束问题的一类直观方法是将其转化为若干个标量的不等式约束。例如，将线性矩阵不等式约束 $F\left(x\right)\succcurlyeq 0$ 转化为 $F\left(x\right)$ 的顺序主子式大于等于0，进而将矩阵型约束转化为标量型约束，继而可以采用传统的标量型松弛方法进行求解，如采用经典的Lasserre松弛方法 [3] 进行求解。Lasserre松弛方法是基于半定规划松弛的，通过构造一系列半定规划松弛问题，可以无限逼近原多项式优化问题的全局最优值。Lasserre证明了当阿基米德条件成立时，该松弛方法具有渐进收敛性。但是，当 $F\left(x\right)$ 的维数很大时，此类方法可能会导致大量的高次约束，导致松弛问题的规模过大，从而使得不适合采用标量型松弛方法进行求解。因此，许多学者致力于直接研究带有多项式矩阵约束的优化问题。例如，Hol等 [4] [5] 和Kojima [6] 对多项式矩阵平方和分解进行了深入的研究。为了求解非凸多项式矩阵不等式优化问题，Henrion等 [7] 利用Hol等 [4] [5] 和Kojima [6] 关于非负多项式矩阵平方和分解的矩理论的研究成果，将标量型松弛方法推广到了求解非凸多项式矩阵不等式优化问题的矩阵型松弛方法。Henrion等 [7] 指出当阿基米德条件成立时，该矩阵型松弛方法具有渐进收敛性。

此外，Nie [8] 给出了利用半定松弛方法来计算一个多项式的局部极小值的层次结构算法。下面我们以在无约束情况下求一个多项式的局部极小值为例对该方法进行说明。具体方法如下：任意给定一个多项式 $f\left(x\right)$ ，利用最优性条件 $\nabla f\left(x\right)=0$ ，构造半定松弛问题并定义H-最小值点以及H-最小值，然后通过构造并求解一系列半定松弛问题，可以将全部的H-最小值点以及H-最小值计算出来，最后可以从这些值中找出 $f\left(x\right)$ 的全体局部极小值，且Nie证明了构造的这一系列松弛问题都具有有限收敛性。

利用线性矩阵不等式的性质 [9] [10]，我们可以将线性矩阵不等式半定互补问题进行等价改写，得到与该问题等价的矩阵形式的问题，进而可以采用矩阵型松弛方法进行求解。进一步，当线性矩阵不等式半定互补问题的实数解个数有限时，利用层次结构算法的思想可以计算出该问题的全部实数解。

符号说明：符号 $\mathbb{N}$ ( $ℝ$ )表示非负整数集(实数集)。符号 ${\mathbb{N}}^n$ 表示n维非负整数向量集合。符号 ${ℝ}^n$ 表示n维实向量集合。对于 $x\in {ℝ}^n$ ， $x_i$ 表示 $x$ 中的第i个元素，即 $x=\left(x_1,x_2,\cdots ,x_n\right)$ 。对于 $\alpha \in {\mathbb{N}}^n$ ，令 $\left|\alpha \right|={\alpha }_1+{\alpha }_2+\cdots +{\alpha }_n$ 。对于 $x\in {ℝ}^n$ 及 $\alpha \in {\mathbb{N}}^n$ ，令 $x^{\alpha }=x_1^{{\alpha }_1}{x}_2^{{\alpha }_2}\cdots x_n^{{\alpha }_n}$ 。对于整数 $n>0$ ， $\left[n\right]$ 表示集合 $\left\{1,2,\cdots ,n\right\}$ 。对于任意整数 $t\in ℝ$ ， $\boxed{t}$ 表示不小于t的最小整数。上标T表示矩阵或向量的转置。对于多项式 $f\left(x\right)$ ， $\text{deg}\left(f\left(x\right)\right)$ 表示它的次数。符号 ${ℝ}^{m\times m}$ 表示 $m\times m$ 维实对称矩阵空间。符号 $ℝ\left[x\right]=ℝ\left[x_1,x_2,\cdots ,x_n\right]$ 表示在实数域上以 $x$ 为自变量的多项式环，即任意 $f\left(x\right)\in ℝ\left[x\right]$ 是一个以 $x$ 为自变量的多项式。令 $ℝ{\left[x\right]}_k$ 是指由次数至多为k的 $ℝ\left[x\right]$ 中的多项式构成的集合，即 $ℝ{\left[x\right]}_k=\left\{f\left(x\right)\in ℝ\left[x\right]:\mathrm{deg}\left(f\left(x\right)\right)\le k\right\}$ 。令 $\mathcal{S}{\left[x\right]}^m$ 表示 $m\times m$ 维实对称多项式矩阵空间，即任意 $F\left(x\right)\in \mathcal{S}{\left[x\right]}^m$ ， $F\left(x\right):{ℝ}^n\mapsto \mathcal{S}{\left[x\right]}^m$ 是一个对称多项式矩阵且每一项 $F_{i,j}\left(x\right)\in ℝ\left[x\right]\left(1\le i,j\le m\right)$ 是一个以 $x$ 为自变量的多项式。如果对于任意非零向量 $u\in {ℝ}^m$ ，二次型 $u^{\text{T}}F\left(x\right)u$ 是非负的，则称矩阵 $F\left(x\right)\in \mathcal{S}{\left[x\right]}^m$ 是正定的，记作 $F\left(x\right)\succcurlyeq 0$ 。一个多项式矩阵 $F\left(x\right)$ 的次数 $\text{deg}\left(F\left(x\right)\right)$ 是矩阵中全部元素所代表的多项式中的最大次数，即 $\text{deg}\left(F\left(x\right)\right)={\text{max}}_{i,j}\text{deg}\left(F_{i,j}\left(x\right)\right)$ 。令 $\langle A,B\rangle =\text{trace}\left(A^{\text{T}}B\right)$ 表示两个向量或矩阵的标准内积。

2. 线性矩阵不等式半定互补问题的等价改写

本节研究如何等价地改写线性矩阵不等式半定互补问题：给定矩阵 $F_1\left(x\right)$ ， $F_2\left(x\right)\in \mathcal{S}{\left[x\right]}^m$ ，寻求解 $x\in {ℝ}^n$ 满足 $F_1\left(x\right)\succcurlyeq 0$ ， $F_2\left(x\right)\succcurlyeq 0$ ， $\langle F_1\left(x\right),F_2\left(x\right)\rangle =0$ ，其中 $F_i\left(x\right)\succcurlyeq 0,i=1,2$ 是线性矩阵不等式。记该问题的全体实数解构成的集合为 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ ，则

$S\left(F_1\left(x\right),F_2\left(x\right)\right)=\left\{x\in {ℝ}^n:F_1\left(x\right)\succcurlyeq 0,F_2\left(x\right)\succcurlyeq 0,\text{trace}\left(F_1\left(x\right)F_2\left(x\right)\right)=0\right\}.$ (1)

根据线性代数的知识可以得到：对于两个半正定矩阵 $A,B$ ， $\text{trace}\left(AB\right)=0$ 当且仅当 $AB=0$ 。因此，线性矩阵不等式半定互补问题中的互补条件等价于矩阵型等式约束，即：

$S\left(F_1\left(x\right),F_2\left(x\right)\right)=\left\{x\in {ℝ}^n:F_1\left(x\right)\succcurlyeq 0,F_2\left(x\right)\succcurlyeq 0,F_1\left(x\right)F_2\left(x\right)=0\right\}.$ (2)

令 $H\left(x\right)=F_1\left(x\right)F_2\left(x\right)$ 是 $m\times m$ 维多项式矩阵，其中 $F_i\left(x\right)={\left[f_{st}^{\left(i\right)}\left(x\right)\right]}_{s,t=1,2,\cdots ,m},i=1,2$ ，则 $H_{ij}\left(x\right)={\displaystyle \underset{k=1}{\overset{m}{\sum }}{f}_{ik}^{\left(1\right)}\left(x\right)f_{kj}^{\left(2\right)}}\left(x\right),1\le i,j\le m$ 。令标量多项式 $h_{\left(i-1\right)m+j}\left(x\right)=H_{ij}\left(x\right),1\le i,j\le m$ ，则有多项式组 $h\left(x\right)=\left(H_{11}\left(x\right),H_{12}\left(x\right),\cdots ,H_{1m}\left(x\right),H_{21}\left(x\right),H_{22}\left(x\right),\cdots ,H_{mm}\left(x\right)\right)=\left(h_1\left(x\right),h_2\left(x\right),\cdots ,h_{m^2}\left(x\right)\right)$ 。因此，可以得到

$S\left(F_1\left(x\right),F_2\left(x\right)\right)=\left\{x\in {ℝ}^n:F_1\left(x\right)\succcurlyeq 0,F_2\left(x\right)\succcurlyeq 0,h_z\left(x\right)=0,z=1,2,\cdots ,m^2\right\}.$ (3)

由公式(3)可以看出，集合 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 的定义中涉及了线性矩阵不等式约束 $F_1\left(x\right)\succcurlyeq 0$ 和 $F_2\left(x\right)\succcurlyeq 0$ 。此前，对于这类约束的优化问题的求解较为困难。事实上，一类处理方法是我们可以将线性矩阵不等式约束转化为标量型约束，从而用一系列标量型不等式约束来描述集合 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 。例如：当 $F\left(x\right)$ 的维数大于2时， $F\left(x\right)\succcurlyeq 0$ 当且仅当它的各阶顺序主子式大于等于零 [11]，即：对于 $i=1,2$ ， $F_i\left(x\right)\succcurlyeq 0$ 等价于 $\left\{g_{ij}\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m\right\}$ ，其中 $g_{ij}\left(x\right)$ 为 $F_i\left(x\right)$ 的第j阶顺序主子式，则可以得到 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 的等价形式：

$S\left(F_1\left(x\right),F_2\left(x\right)\right)=\left\{x\in {ℝ}^n:h_z\left(x\right)=0,g_{ij}\left(x\right)\ge 0,z=1,2,\cdots ,m^2;i=1,2;j=1,2,\cdots ,m\right\}.$ (4)

特别地，当 $F\left(x\right)$ 的维数等于2时， $F\left(x\right)\succcurlyeq 0$ 当且仅当 $\text{trace}\left(F\left(x\right)\right)\ge 0$ 且 $\text{det}\left(F\left(x\right)\right)\ge 0$ ，则

$\begin{array}{c}S\left(F_1\left(x\right),F_2\left(x\right)\right)=\{x\in {ℝ}^n:\text{trace}\left(F_1\left(x\right)\right)\ge 0,\det \left(F_1\left(x\right)\right)\ge 0,\text{trace}\left(F_2\left(x\right)\right)\ge 0,\\ \det \left(F_2\left(x\right)\right)\ge 0,h_z\left(x\right)=0,z=1,2,\cdots ,m^2\}.\end{array}$

由公式(4)中可以看出，当原问题的规模较大时，标量化改写方法会引入大量的高次约束，导致松弛问题的规模过大，从而变得难以求解。因此我们考虑保留线性矩阵不等式半定互补问题的矩阵型约束进行计算。

由文献 [9] 可知，线性矩阵不等式约束序列 $F_1\left(x\right)\succcurlyeq 0,F_2\left(x\right)\succcurlyeq 0,\cdots ,F_k\left(x\right)\succcurlyeq 0$ 等价于一个线性矩阵不等式约束 $\text{diag}\left(F_1\left(x\right),F_2\left(x\right),\cdots ,F_k\left(x\right)\right)\succcurlyeq 0$ 。

令 $Q\left(x\right)=\text{diag}\left(F_1\left(x\right),F_2\left(x\right)\right)$ ，则

(5)

由公式(4)和(5)可知，通过随机地选取一个恰当的多项式 $f\left(x\right)$ 作为目标函数，我们可以将原线性矩阵不等式半定互补问题的解集分别转化为下列两种形式的多项式优化问题的可行域，其中标量型多项式优化问题为：

$\begin{array}{l}\underset{x\in {ℝ}^n}{\min }f\left(x\right)\\ \text{s}\text{.t}\text{.}{g}_{ij}\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m,\\ h_z\left(x\right)=0,z=1,2,\cdots ,m^2,\end{array}$ (6)

矩阵型多项式优化问题为：

(7)

显然， $x\in S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 当且仅当 $x$ 是问题(6)或(7)的一个可行解。

3. 求解全体实数解的层次结构方法

本节中，我们将研究线性矩阵不等式半定互补问题的层次结构方法，假设实数解的个数有限，我们将给出求解该问题的全体实数解的方法。

首先，我们考虑在问题(6)和问题(7)中选取目标函数 $f\left(x\right)$ 。事实上，通过适当地选取，多项式 $f\left(x\right)$ 在 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 的不同的解处取得不同的函数值，将这些值依次排列 $f_1<f_2<\cdots <f_N$ ，称 $f_i$ 是第i个可行值。

令第i个解集合为 $S_i=S\left(F_1\left(x\right),F_2\left(x\right)\right)\cap \left\{x\in {ℝ}^n:f\left(x\right)=f_i\right\}$ ，则解集合 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 能被划分为有限多个互不相交的非空子集，即 $S\left(F_1\left(x\right),F_2\left(x\right)\right)=S_1\cup S_2\cup \cdots \cup S_N$ 。

若我们可以将解集合 $S_i,i=1,2,\cdots ,N$ 全部求解出来，则可以计算出原问题的全部实数解。下面我们分别构造与每一个子集合 $S_i$ 相对应的优化子问题，进而求解出所有的 $S_i$ 。此处类比Nie在文献 [8] 中求解多项式的局部极小值的层次结构算法完成这项工作。

显然，第一个解集合 $S_1$ 是第2节中构造的两个多项式优化问题(6)和(7)的最优解的集合，即，标量型优化问题

$\begin{array}{l}{f}_1^{1*}:=\underset{x\in {ℝ}^n}{\min }f\left(x\right)\\ \text{s}\text{.t}\text{.}{g}_{ij}\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m,\\ h_z\left(x\right)=0,z=1,2,\cdots ,m^2,\end{array}$ (8)

和等价的矩阵型优化问题

(9)

如果问题(8)或问题(9)是可解的，则可以得到问题的最优值，在不需要对比两种方法时，为了方便叙述都将其记为 $f_1$ ，称问题(8)或问题(9)是第一个子问题。

获得了可行值 $f_1$ 和相应的可行集 $S_1$ 之后，进一步构造并求解下一个子问题以求解第二个可行集 $S_2$ 。依次进行下去，每一步在求解出可行值 $f_t$ 和相应的可行集 $S_t$ 之后，构造并求解下一个子问题以求解可行集 $S_{t+1}$ ，最终可以求解出全部可行集 $S_1,S_2,\cdots ,S_N$ 。具体的方法如下：在获得了前t个可行值 $f_1,f_2,\cdots ,f_t$ 和相应的可行集 $S_1,S_2,\cdots ,S_t$ 之后，构造第 $\left(t+1\right)$ 个子问题，其中第 $\left(t+1\right)$ 个标量型子问题为：

$\begin{array}{l}{f}_{t+1}^{1*}=\underset{x\in {ℝ}^n}{\min }f\left(x\right)\\ \text{s}\text{.t}\text{.}{g}_{ij}\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m,\\ h_z\left(x\right)=0,z=1,2,\cdots ,m^2,\\ f\left(x\right)\ge f_t+\epsilon ,\end{array}$ (10)

第 $\left(t+1\right)$ 个矩阵型子问题为：

(11)

其中问题(10)和问题(11)中使用的参数 $ϵ$ 是满足

$0<ϵ<f_{t+1}-f_t,$ (12)

的一个正数。当第 $\left(t+1\right)$ 个子问题可解时，无论标量型问题(10)还是等价的矩阵型问题(11)，在不需要对比两种方法时，为了方便叙述都将其最优值记为 $f^{>}\left(f_t+ϵ\right)$ ，在需要对比两种方法时，我们将这两个最优值分别记为 $f_{t+1}^{1*}$ 和 $f_{t+1}^{2*}$ 。

引理3.1 $f^{>}\left(f_t+ϵ\right)=f_{t+1}$ 当且仅当参数 $ϵ$ 满足条件(12)。

证明：如果 $ϵ$ 满足条件(12)，则 $f_{t+1}>f_t+\epsilon$ 。又因为 $f^{>}\left(f_t+ϵ\right)$ 是 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 上大于或等于 $f_t+\epsilon$ 的最小可行值，则 $f^{>}\left(f_t+ϵ\right)=f_{t+1}$ 。反之，如果 $f^{>}\left(f_t+ϵ\right)=f_{t+1}$ ，则表明在 $\left(f_t+\epsilon ,+\infty \right)$ 上存在可行值，所以，若 $f_{t+1}$ 存在，则必有 $f_{t+1}>f_t+\epsilon$ ，即参数 $ϵ$ 满足条件(12)。

但是，考虑到 $f_{t+1}$ 是否存在是未知的，即便知道其存在，也很难提前去确定它的值。在实际计算过程中，参数 $ϵ$ 的取值不宜过小，以免在数值计算过程中导致结果产生较大偏差。因此，我们需要考虑如何找到合适的参数 $ϵ$ 来构造第 $\left(t+1\right)$ 个子问题(10)或(11)。为了解决这一问题，我们分别构造第 $\left(t+1\right)$ 个判定优化子问题：即标量型判定问题

$\begin{array}{l}\underset{x\in {ℝ}^n}{\max }f\left(x\right)\\ \text{s}\text{.t}\text{.}{g}_{ij}\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m,\\ h_z\left(x\right)=0,z=1,2,\cdots ,m^2,\\ f\left(x\right)\le f_t+\epsilon ,\end{array}$ (13)

和矩阵型判定问题

(14)

明显地，上述两类判定优化子问题都是可解的，因为任意 $x\in S_1\cup S_2\cup \cdots \cup S_t$ 都是它的可行解，我们将问题(13)和问题(14)的最优值记为 $f^{<}\left(f_t+ϵ\right)$ ，则有下列结果成立。

引理3.2对于 $ϵ>0$ ， $f^{<}\left(f_t+ϵ\right)=f_t$ 当且仅当 $ϵ<f_{t+1}-f_t$ 。

证明：如果 $ϵ<f_{t+1}-f_t$ ，则 $f_{t+1}>f_t+\epsilon$ ，由问题(13)或(14)可知， $f^{<}\left(f_t+ϵ\right)$ 是 $S\left(F_1\left(x\right),F_2\left(x\right)\right)$ 上小于或等于 $f_t+\epsilon$ 的最大可行值，则 $f^{<}\left(f_t+ϵ\right)=f_t$ 。反之，如果 $f^{<}\left(f_t+ϵ\right)=f_t$ ，则表明在 $\left(f_t,f_t+\epsilon \right)$ 上不存在可行值，所以，如果 $f_{t+1}$ 存在，则必有 $f_{t+1}>f_t+\epsilon$ ，即 $ϵ<f_{t+1}-f_t$ 成立。

事实上， $ϵ$ 的值不能取太小，例如0.05是符合计算上的要求的。在问题(13)或问题(14)中，令 $ϵ=0.05$ 。如果 $f^{<}\left(f_t+ϵ\right)>f_t$ ，我们将 $ϵ$ 减小至 $ϵ/5$ ，再次求解问题(13)或问题(14)，重复这个过程直至 $f^{<}\left(f_t+ϵ\right)=f_t$ ，我们即可得到合适的参数 $ϵ$ 的值。

利用引理3.2的结果，我们可以确定适当的参数 $ϵ$ ，利用该参数 $ϵ$ 和当前可行值 $f_t$ ，我们可以构造下一个子问题，参考公式(10)和(11)。如果下一个子问题不可行，则说明下一个可行值 $f_{t+1}$ 不存在，即 $f_t$ 是线性矩阵不等式半定互补问题的所有的可行值中的最大值。否则，若 $f_{t+1}$ 存在且 $ϵ$ 满足条件(12)，通过求解第 $\left(t+1\right)$ 个子问题可以得到下一个可行集 $S_{t+1}$ 。至此，我们给出了线性矩阵不等式半定互补问题的层次结构方法。若我们能够求解出标量型子问题序列(10)和矩阵型多项式优化子问题序列(11)，以及它们对应的判定问题序列(13)和(14)，即可求解出原问题的全部实数解。

4. 半定松弛算法及其收敛性

在前两节的内容中，我们利用层次结构方法将线性矩阵不等式半定互补问题转化为标量型多项式优化子问题序列(8)，(10)，(13)或矩阵型多项式优化子问题序列(9)，(11)，(14)。本节中，我们对这两类多项式优化子问题分别给出半定松弛求解方法，从而完成对原问题的求解。

4.1. 求解标量型多项式优化子问题的半定松弛方法

对于标量型子问题序列(8)，(10)和(13)，我们采用Lasserre在文献 [3] 中提出的经典半定松弛方法进行求解。

首先，我们需要定义局部化矩矩阵。令 $v\left(x\right)={\left[1,x_1,\cdots ,x_n,x_1^2,x_1{x}_2,\cdots ,x_1^k,\cdots ,x_n^k,\cdots \right]}^{\text{T}}$ 和 $v_k\left(x\right)={\left[1,x_1,\cdots ,x_n,x_1^2,x_1{x}_2,\cdots ,x_1^k,\cdots ,x_n^k\right]}^{\text{T}}$ 分别是 $ℝ\left[x\right]$ 和 $ℝ{\left[x\right]}_k$ 的一组基，其中 $v_k\left(x\right)\in ℝ{\left[x\right]}^{s\left(k\right)}$ ， $s\left(k\right)=\left(\begin{array}{c}n+k\\ n\end{array}\right)$ 是它的维数。令 $y={\left\{y_{\alpha }\right\}}_{\alpha \in {\mathbb{N}}^n}\in {ℝ}^{{\mathbb{N}}_d^n}$ 是一个矩序列， $p\left(x\right)$ 是任意一个多项式，定义Riesz函数 ${\mathcal{F}}_y:ℝ\left[x\right]\mapsto ℝ{\left[x\right]}_k$ ： ${\mathcal{F}}_y\left(p\right)={\displaystyle \underset{\alpha \in {\mathbb{N}}^n}{\sum }{p}_{\alpha }}{y}_{\alpha }$ 。如果给定一个截断矩序列 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ ，定义关于 $y$ 的k阶矩矩阵为 $M_k\left(y\right)={\mathcal{F}}_y\left(v_k\left(x\right)v_k{\left(x\right)}^{\text{T}}\right)$ 。对于标量多项式 $l_t\left(x\right)\in ℝ\left[x\right]$ ，定义由截断矩序列 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ 生成的关于多项式 $l_t\left(x\right)\in ℝ\left[x\right]$ 的k阶局部化矩矩阵为 $L_{l_t\left(x\right)}^{\left(k\right)}\left(y\right)={\mathcal{F}}_y\left(l_t\left(x\right)v_{d_1}\left(x\right)v_{d_1}{\left(x\right)}^{\text{T}}\right)$ ，其中 $d_1=k-\boxed{\mathrm{deg}\left(l_t\right)/2}$ 。

如果一个截断矩序列 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ 满足下面的条件：

$\begin{array}{l}{M}_k\left(y\right)\succcurlyeq 0,L_{g_{ij}\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,L_{h_z\left(x\right)}^{\left(k\right)}\left(y\right)=0,L_{l_t\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,\\ s=\text{rank}{M}_k\left(y\right)=\text{rank}{M}_{k-d}\left(y\right),\end{array}$ (15)

其中 $d=\text{max}\left\{1,\boxed{\text{deg}\left(g_{ij}\right)/2},\boxed{\text{deg}\left(h_z\right)/2},\boxed{\text{deg}\left(l_t\right)/2}\right\}$ ，则称 $y$ 满足平滑条件。若 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ 满足平滑条件(15)，则 $y$ 在半代数集

$\left\{x\in {ℝ}^n:g_{ij}\left(x\right)\ge 0,h_z\left(x\right)=0,l_t\left(x\right)\ge 0,i=1,2;j=1,2,\cdots ,m;z=1,2,\cdots ,m^2\right\}$

上存在唯一的一个s-原子表示测度。

接下来，我们给出矩阵型子问题序列(8)，(10)和(13)的求解方法。例如，采用Lasserre松弛方法求解问题(10)时，创建问题(10)的k阶松弛如下：

$\begin{array}{l}{r}_{t+1,k}^{1*}=\min L_f\left(y\right)\\ \text{s}\text{.t}\text{.}{M}_k\left(y\right)\succcurlyeq 0,\\ L_{g_{ij}\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,i=1,2;j=1,2,\cdots ,m,\\ L_{h_z\left(x\right)}^{\left(k\right)}\left(y\right)=0,z=1,2,\cdots ,m^2,\\ L_{l_t\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,\\ y_0=1,\\ y\in {ℝ}^{{\mathbb{N}}_{2k}^n},\end{array}$ (16)

其中 $l_t\left(x\right)=f\left(x\right)-f_t-ϵ$ ， $k=k_0,k_0+1,\cdots$ ，且

$k_0=\text{max}\left\{1,\boxed{\text{deg}\left(f\right)/2},\boxed{\text{deg}\left(g_{ij}\right)/2},\boxed{\text{deg}\left(h_z\right)/2},\boxed{\text{deg}\left(l_t\right)/2}\right\}$

是最低松弛阶数。

令 $I_{2k}\left(h\right)$ 是由多项式组 $h\left(x\right)$ 生成的理想的第2k阶截断，即

$I_{2k}\left(h\right)=h_1\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_1\right)}+h_2\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_2\right)}+\cdots +h_{m^2}\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_{m^2}\right)},$

$Q_k\left(g_{ij},l_t\right)$ 是由 $g_{ij}\left(x\right)$ 和 $l_t\left(x\right)$ 生成的二次模块的第k阶截断，即

$Q_k\left(g_{ij},l_t\right)={{\displaystyle \sum \left[x\right]}}_{2k}+g_{ij}\left(x\right)\cdot {{\displaystyle \sum \left[x\right]}}_{2k-\mathrm{deg}\left(g_{ij}\right)}+l_t\left(x\right)\cdot {{\displaystyle \sum \left[x\right]}}_{2k-\mathrm{deg}\left(l_t\right)}.$

问题(16)的对偶是：

$\begin{array}{l}{\lambda }_{t+1,k}^{1*}=\underset{x\in {ℝ}^n}{\max }\gamma \\ \text{s}\text{.t}.f\left(x\right)-\gamma \in I_{2k}\left(h\right)+Q_k\left(g_{ij},l_t\right).\end{array}$ (17)

半定规划问题(16)和问题(17)是对偶问题。因此，对于多项式优化问题(10)，可以采用原对偶内点方法 [12] 求解相应的松弛问题。通过求解一系列松弛问题(16)，可以得到序列 $\left\{r_{t+1,k}^{1*}\right\}$ ，随着松弛阶数增大， $\left\{r_{t+1,k}^{1*}\right\}$ 逐渐逼近于多项式优化问题(10)的最优值。问题(8)和(13)可以采用类似的方法进行求解。

4.2. 求解矩阵型多项式优化子问题的半定松弛方法

对于矩阵型子问题序列(9)，(11)和(14)，我们采用矩阵型半定松弛方法进行求解。类似于标量型Lasserre半定松弛方法，我们参考Henrion等在文献 [7] 中提出的矩阵型松弛方法。首先，定义由矩阵 $Q\left(x\right)$ 生成的局部化矩矩阵为 $L_{Q\left(x\right)}^{\left(k\right)}\left(y\right)={\mathcal{F}}_y\left(v_{d_2}\left(x\right)v_{d_2}{\left(x\right)}^{\text{T}}\otimes Q\left(x\right)\right)$ ，其中 $d_2=k-\boxed{\text{deg}\left(Q\right)/2}$ ， $\otimes$ 表示Kronecker积。

如果一个截断矩序列 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ 满足下面的条件：

$\begin{array}{l}{M}_k\left(y\right)\succcurlyeq 0,L_{Q\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,L_{h_z\left(x\right)}^{\left(k\right)}\left(y\right)=0,L_{l_t\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,\\ s=\text{rank}{M}_k\left(y\right)=\text{rank}{M}_{k-d}\left(y\right),\end{array}$ (18)

其中 $d=\text{max}\left\{1,\boxed{\text{deg}\left(Q\right)/2},\boxed{\text{deg}\left(h_z\right)/2},\boxed{\text{deg}\left(l_t\right)/2}\right\}$ ，则称 $y$ 满足平滑条件。若 $y\in {ℝ}^{{\mathbb{N}}_{2k}^n}$ 满足平滑条件(18)，则 $y$ 在半代数集

$\left\{x\in {ℝ}^n:Q\left(x\right)\succcurlyeq 0,h_z\left(x\right)=0,l_t\left(x\right)\ge 0,z=1,2,\cdots ,m^2\right\}$

上存在唯一的一个s-原子表示测度。

接下来，我们给出矩阵型子问题序列(9)，(11)和(14)的求解方法。例如，采用矩阵型松弛方法求解问题(11)时，创建问题(11)的k阶松弛如下：

$\begin{array}{l}{r}_{t+1,k}^{2*}=\min L_f\left(y\right)\\ \text{s}\text{.t}.M_k\left(y\right)\succcurlyeq 0,\\ L_{Q\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,\\ L_{h_z\left(x\right)}^{\left(k\right)}\left(y\right)=0,z=1,2,\cdots ,m^2,\\ L_{l_t\left(x\right)}^{\left(k\right)}\left(y\right)\succcurlyeq 0,\\ y_0=1,\\ y\in {ℝ}^{{\mathbb{N}}_{2k}^n},\end{array}$ (19)

其中 $l_t\left(x\right)=f\left(x\right)-f_t-ϵ$ ， $k=k_0,k_0+1,\cdots$ ，且

$k_0=\text{max}\left\{1,\boxed{\text{deg}\left(f\right)/2},\boxed{\text{deg}\left(Q\right)/2},\boxed{\text{deg}\left(h_z\right)/2},\boxed{\text{deg}\left(l_t\right)/2}\right\}$

是最低松弛阶数。

令 $I_{2k}\left(h\right)$ 是由多项式组 $h\left(x\right)$ 生成的理想的第2k阶截断，即

$I_{2k}\left(h\right)=h_1\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_1\right)}+h_2\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_2\right)}+\cdots +h_{m^2}\cdot ℝ{\left[x\right]}_{2k-\mathrm{deg}\left(h_{m^2}\right)},$

$Q_k\left(Q,l_t\right)$ 是由 $Q\left(x\right)$ 和 $l_t\left(x\right)$ 生成的二次模块的第k阶截断，即

$\begin{array}{l}{Q}_k\left(Q,l_t\right)={{\displaystyle \sum \left[x\right]}}_{2k}+l_t\left(x\right)\cdot {{\displaystyle \sum \left[x\right]}}_{2k-\mathrm{deg}\left(l_t\right)}\\ +\left\{\langle R,Q\rangle :R\in {{\displaystyle \sum \left[x\right]}}^{m\times m},\mathrm{deg}\left(R_{\left(i,j\right)}{Q}_{\left(i,j\right)}\right)\le 2k,\forall i,j\in \left[m\right]\right\}.\end{array}$

问题(19)的对偶是：

$\begin{array}{l}{\lambda }_{t+1,k}^{2*}=\underset{x\in {ℝ}^n}{\max }\gamma \\ \text{s}\text{.t}\text{.}f\left(x\right)-\gamma \in I_{2k}\left(h\right)+Q_k\left(Q,l_t\right).\end{array}$ (20)

半定规划问题(19)和问题(20)是对偶问题。因此，对于多项式优化问题(11)，可以采用原对偶内点方法 [12] 求解相应的松弛问题。通过求解一系列松弛问题(19)，可以得到序列 $\left\{r_{t+1,k}^{2*}\right\}$ ，随着松弛阶数增大， $\left\{r_{t+1,k}^{2*}\right\}$ 逐渐逼近于多项式优化问题(11)的最优值。问题(9)和(14)可以采用类似的方法进行求解。

4.3. 算法

下面，我们给出求解线性矩阵不等式半定互补问题的全体实数解的半定松弛算法。

算法4.3.1：求解线性矩阵不等式半定互补问题的标量型算法

步骤1. 随机地选取一个多项式 $f\left(x\right)$ 作为目标多项式。令 $t:=0$ ， $k:=k_0$ 。符号 $S\left(F_1,F_2\right)$ 表示线性矩阵不等式半定互补问题的解集， $\mathcal{F}$ 表示可行值构成的集合。

步骤2. 若问题(8)的松弛问题是不可行的，则令 $S\left(F_1,F_2\right)=\varnothing$ ，停止。否则，求解问题(8)的松弛问题，得到最优解 $y_{1*}^{1,k}$ 和最优值 $r_{1,k}^{1*}$ 。

步骤3. 若存在某个 $s\in \left[k_0,k\right]$ ，截断矩序列 $y_{1*}^{1,k}{|}_{2s}$ 满足平滑条件(15)，则令 $S\left(F_1,F_2\right):=S_1$ ，其中 $S_1$ 是问题(8)的最优解的集合。令 $\mathcal{F}=\left\{r_{1,k}^{1*}\right\}$ ， $k:=k_0$ ， $t:=t+1$ ，转向步骤4。若不存在这样的s，则令 $k:=k+1$ ，转向步骤2。

步骤4. 令 $ϵ=0.05$ ，计算问题(13)的最优值 $f^{<}\left(f_t+ϵ\right)$ 。若 $f^{<}\left(f_t+ϵ\right)>f_t$ ，则令 $ϵ=ϵ/5$ ，再次求解问题(13)，重复这个过程直至 $f^{<}\left(f_t+ϵ\right)=f_t$ 。

步骤5. 求解问题(10)的松弛问题。若它不可行，则没有更多的实数解，令 $S\left(F_1,F_2\right)=S$ ，其中 $S=S_1\cup S_2\cup \cdots \cup S_t$ ，停止；否则，求解问题(10)的松弛问题得到最优解 $y_{1*}^{t+1,k}$ 和最优值 $r_{t+1,k}^{1*}$ 。

步骤6. 若存在某个 $s\in \left[k_0,k\right]$ ，截断矩序列 $y_{1*}^{t+1,k}{|}_{2s}$ 满足平滑条件(15)，则更新 $S:=S\cup S_{t+1}$ ，其中 $S_{t+1}$ 是问题(10)的最小值点的集合。令 $\mathcal{F}:=\mathcal{F}\cup \left\{r_{t+1,k}^{1*}\right\}$ ， $k:=k_0$ ， $t:=t+1$ ，转向步骤4。若不存在这样的s，则令 $k:=k+1$ ，转向步骤5。

算法4.3.2：求解线性矩阵不等式半定互补问题的矩阵型算法

步骤1. 随机地选取一个多项式 $f\left(x\right)$ 作为目标多项式。令 $t:=0$ ， $k:=k_0$ 。符号 $S\left(F_1,F_2\right)$ 表示线性矩阵不等式半定互补问题的解集， $\mathcal{F}$ 表示可行值构成的集合。

步骤2. 若问题(9)的松弛问题是不可行的，则令 $S\left(F_1,F_2\right)=\varnothing$ ，停止。否则，求解问题(9)的松弛问题，得到最优解 $y_{2*}^{1,k}$ 和最优值 $r_{1,k}^{2*}$ 。

步骤3. 若存在某个 $s\in \left[k_0,k\right]$ ，截断矩序列 $y_{2*}^{1,k}{|}_{2s}$ 满足平滑条件(18)，则令 $S\left(F_1,F_2\right):=S_1$ ，其中 $S_1$ 是问题(9)的最优解的集合。令 $\mathcal{F}=\left\{r_{1,k}^{2*}\right\}$ ， $k:=k_0$ ， $t:=t+1$ ，转向步骤4。若不存在这样的s，则令 $k:=k+1$ ，转向步骤2。

步骤4. 令 $ϵ=0.05$ ，计算问题(14)的最优值 $f^{<}\left(f_t+ϵ\right)$ 。若 $f^{<}\left(f_t+ϵ\right)>f_t$ ，则令 $ϵ=ϵ/5$ ，再次求解问题(14)，重复这个过程直至 $f^{<}\left(f_t+ϵ\right)=f_t$ 。

步骤5.求解问题(11)的松弛问题。若它不可行，则没有更多的实数解，令 $S\left(F_1,F_2\right)=S$ ，其中 $S=S_1\cup S_2\cup \cdots \cup S_t$ ，停止；否则，求解问题(11)的松弛问题得到最优解 $y_{2*}^{t+1,k}$ 和最优值 $r_{t+1,k}^{2*}$ 。

步骤6.若存在某个 $s\in \left[k_0,k\right]$ ，截断矩序列 $y_{2*}^{t+1,k}{|}_{2s}$ 满足平滑条件(18)，则更新 $S:=S\cup S_{t+1}$ ，其中 $S_{t+1}$ 是问题(11)的最小值点的集合。令 $\mathcal{F}:=\mathcal{F}\cup \left\{r_{t+1,k}^{2*}\right\}$ ， $k:=k_0$ ， $t:=t+1$ ，转向步骤4。若不存在这样的s，则令 $k:=k+1$ ，转向步骤5。

对于Lasserre松弛方法 [3]，Lasserre证明了在阿基米德条件成立时，该松弛方法具有渐进收敛性。Henrion等 [7] 利用非负多项式矩阵平方和分解的相关理论结果，将标量型松弛方法进行推广，得到了矩阵型松弛方法，并证明了在阿基米德条件成立时，该矩阵型松弛方法同样具有渐进收敛性。因此，通过文献 [3] 和 [7] 的相关结论，可以得到算法4.3.1和算法4.3.2的收敛性定理如下：

定理4.3.3假设线性矩阵不等式半定互补问题的实数解个数有限，则当松弛阶数 $k\to \infty$ 时， $r_{t,k}^{1*}\to f_t^{1*},{\lambda }_{t,k}^{1*}\to f_t^{1*},r_{t,k}^{2*}\to f_t^{2*},{\lambda }_{t,k}^{2*}\to f_t^{2*}$ 。

根据定理4.3.3可知：当线性矩阵不等式半定互补问题的实数解个数有限时，算法4.3.1和4.3.2具有渐进收敛性。

5. 数值实验

本节将给出求解线性矩阵不等式半定互补问题的数值实验结果。下表中的运行时间是每个维度20次的平均运行时间，时间结果保留四位小数。我们分别用下面的两种方法求解线性矩阵不等式半定互补问题：

方法1 (标量型方法)：采用直观的标量化方法，将原问题中的半定约束转化为由顺序主子式组成的约束。将原问题等价转化为子问题序列(8)、(10)和(13)，再采用上一节的算法4.3.1进行求解。

方法2 (矩阵型方法)：保留线性矩阵不等式半定互补问题的矩阵形式进行计算。将原问题等价转化为子问题序列(9)、(11)和(14)，再采用上一节的算法4.3.2进行求解。

在数值实验部分，我们使用的是操作系统是Windows10，内存是8.0 GB RAM，处理器是Inter(R)Core(TM)i58265U@1.80 GHz的笔记本电脑。我们调用Matlab R2014b软件，使用YALMIP [13]，SeDuMi [14] 和Gloptipoly 3 [15] 求解线性矩阵不等式半定互补问题。

例1. 给定多项式矩阵

$F_1\left(x\right)=\left(\begin{array}{cc}{f}_{11}^{\left(1\right)}\left(x\right)& f_{12}^{\left(1\right)}\left(x\right)\\ f_{21}^{\left(1\right)}\left(x\right)& f_{22}^{\left(1\right)}\left(x\right)\end{array}\right)=\left(\begin{array}{cc}1-4x_1& x_1\\ x_1& 4-x_1-x_2\end{array}\right),$

$F_2\left(x\right)=\left(\begin{array}{cc}{f}_{11}^{\left(2\right)}\left(x\right)& f_{12}^{\left(2\right)}\left(x\right)\\ f_{21}^{\left(2\right)}\left(x\right)& f_{22}^{\left(2\right)}\left(x\right)\end{array}\right)=\left(\begin{array}{cc}4-x_2& x_1\\ x_1& x_2\end{array}\right).$

考虑线性矩阵不等式半定互补问题： $F_1\left(x\right)\succcurlyeq 0$ ， $F_2\left(x\right)\succcurlyeq 0$ ， $\langle F_1\left(x\right),F_2\left(x\right)\rangle =0$ 。设目标函数为 $f\left(x\right)=x_2$ 。

记：P_i：求解第i个子问题至达到收敛所需要的松弛次数；S_i：第i个子问题的松弛问题中半定不等式约束的规模。下表中的分号是用于分隔在不同松弛阶数的松弛问题中的数据。为了方便表示，下表中松弛问题的规模S_i用数字表示，例如数字5表示的是在松弛问题中有一个5 × 5的矩阵。我们分别用方法1和方法2求解该问题，求解结果见表1和表2：

表1. 问题的解、计算时间以及第i个子问题的松弛次数

表2. 第i个子问题的松弛问题中半定不等式约束的规模

例2. 给定多项式矩阵

$F_1\left(x\right)=\left(\begin{array}{cc}{f}_{11}^{\left(1\right)}\left(x\right)& f_{12}^{\left(1\right)}\left(x\right)\\ f_{21}^{\left(1\right)}\left(x\right)& f_{22}^{\left(1\right)}\left(x\right)\end{array}\right)=\left(\begin{array}{cc}1-x_1& x_1\\ x_1& 1-x_2\end{array}\right),$

$F_2\left(x\right)=\left(\begin{array}{cc}{f}_{11}^{\left(2\right)}\left(x\right)& f_{12}^{\left(2\right)}\left(x\right)\\ f_{21}^{\left(2\right)}\left(x\right)& f_{22}^{\left(2\right)}\left(x\right)\end{array}\right)=\left(\begin{array}{cc}1-x_2& x_1\\ x_1& 0\end{array}\right).$

考虑线性矩阵不等式半定互补问题： $F_1\left(x\right)\succcurlyeq 0$ ， $F_2\left(x\right)\succcurlyeq 0$ ， $\langle F_1\left(x\right),F_2\left(x\right)\rangle =0$ 。设目标函数为 $f\left(x\right)=x_2$ 。下表中的各个量以及符号所表示的含义与例1相同。我们分别用方法1和方法2求解该问题，求解结果见表3和表4：

表3. 问题的解、计算时间以及第i个子问题的松弛次数

表4. 第i个子问题的松弛问题中半定不等式约束的规模

例3. 给定多项式矩阵

$F_1\left(x\right)=\left(\begin{array}{cccc}{f}_{11}^{\left(1\right)}\left(x\right)& f_{12}^{\left(1\right)}\left(x\right)& f_{13}^{\left(1\right)}\left(x\right)& f_{14}^{\left(1\right)}\left(x\right)\\ f_{21}^{\left(1\right)}\left(x\right)& f_{22}^{\left(1\right)}\left(x\right)& f_{23}^{\left(1\right)}\left(x\right)& f_{24}^{\left(1\right)}\left(x\right)\\ f_{31}^{\left(1\right)}\left(x\right)& f_{32}^{\left(1\right)}\left(x\right)& f_{33}^{\left(1\right)}\left(x\right)& f_{34}^{\left(1\right)}\left(x\right)\\ f_{41}^{\left(1\right)}\left(x\right)& f_{42}^{\left(1\right)}\left(x\right)& f_{43}^{\left(1\right)}\left(x\right)& f_{44}^{\left(1\right)}\left(x\right)\end{array}\right)=\left(\begin{array}{cccc}1-4x_1& x_1& x_1& 0\\ x_1& 0& 4x_1& 0\\ x_1& 4x_1& 1+x_1& x_1\\ 0& 0& x_1& 5x_1\end{array}\right),$