矩陣流形上的數學優化及其應用

(1)

科技部補助專題研究計畫成果報告

期末報告

矩陣流形上的數學優化及其應用(第1年)

計畫類別：個別型計畫

計畫編號： MOST 107-2115-M-006-007-MY2

執行期間： 107年08月01日至108年07月31日

執行單位：國立成功大學數學系暨應用數學所

計畫主持人：林敏雄

計畫參與人員：此計畫無其他參與人員

報告附件：移地研究心得報告

出席國際學術會議心得報告

中　華　民　國　108　年　06　月　24　日

(2)

是Weyl-Horn條件，它表示任意矩陣的特徵值和奇異值之間的關係。

Weyl-Horn的這個結果導致了一個有趣的反問題，即如何構造具所需

特徵值和奇異值的矩陣。在這項工作中，我們提出了來自微分幾何

和不精確牛頓方法的折衷混合技術，用於求解反特徵值和奇異值問

題以及其他具特別特徵的反問題，例如反求之矩陣元素需非負，對

角線元素需固定，甚至預先固定某特定元素。我們從理論上證明了

我們的方法具全局與二次收斂性，這些結果可以從我們提供的數值

例子來看出我們提出的方法的完善性和準確性。在該工作上，鑑於

理論上的興趣，在附錄中我們提供了一個 2乘2實矩陣(甚至是非負

矩陣)具有規定的特徵值，奇異值和主對角線條件的解存在之必要且

充分條件。

中文關鍵詞：反特徵值和奇異值問題，非負矩陣，黎曼不精確牛頓法。

英文摘要： Inverse eigenvalue and singular value problems have been

widely discussed for decades. The well-known result is the

Weyl-Horn condition, which presents the relations between

the eigenvalues and singular values of an arbitrary

matrix. This result by Weyl-Horn then leads to an

interesting inverse problem, i.e., how to construct a

matrix with desired eigenvalues and singular values. In

this work, we do that and more. We propose an eclectic mix

of techniques from differential geometry and the inexact

Newton method for solving inverse eigenvalue and singular

value problems as well as additional desired

characteristics such as nonnegative entries, prescribed

diagonal entries, and even predetermined entries. We show

theoretically that our method converges globally and

quadratically, and we provide numerical examples to

demonstrate the robustness and accuracy of our proposed

method. Having theoretical interest, we provide in the

appendix a necessary and sufficient condition for the

existence of a 2-by-2 real matrix, or even a nonnegative

matrix, with prescribed eigenvalues, singular values, and

main diagonal entries.

英文關鍵詞： Inverse eigenvalue and singular value problems, nonnegative

matrices,

(3)

(will be inserted by the editor)

Riemannian Inexact Newton Method for Structured Inverse

Eigenvalue and Singular Value Problems

Chun-Yueh Chiang · Matthew M. Lin · Xiao-Qing Jin

Received: date / Accepted: date

Abstract Inverse eigenvalue and singular value problems have been widely dis-cussed for decades. The well-known result is the Weyl-Horn condition, which presents the relations between the eigenvalues and singular values of an arbitrary matrix. This result by Weyl-Horn then leads to an interesting inverse problem, i.e., how to con-struct a matrix with desired eigenvalues and singular values. In this work, we do that and more. We propose an eclectic mix of techniques from differential geome-try and the inexact Newton method for solving inverse eigenvalue and singular value problems as well as additional desired characteristics such as nonnegative entries, prescribed diagonal entries, and even predetermined entries. We show theoretically that our method converges globally and quadratically, and we provide numerical ex-amples to demonstrate the robustness and accuracy of our proposed method. Having theoretical interest, we provide in the appendix a necessary and sufficient condition for the existence of a 2 × 2 real matrix, or even a nonnegative matrix, with prescribed eigenvalues, singular values, and main diagonal entries.

Mathematics Subject Classification (2000) 15A29 · 65H17.

Keywords Inverse eigenvalue and singular value problems, nonnegative matrices, Riemannian inexact Newton method.

Center for General Education, National Formosa University, Huwei 632, Taiwan (chiang@nfu.edu.tw). This research was supported in part by the Ministry of Science and Technology of Taiwan under grant 105-2115-M-150-001.

Corresponding author. Department of Mathematics, National Cheng Kung University, Tainan 701, Tai-wan (mhlin@mail.ncku.edu.tw). This research was supported in part by the Ministry of Science and Technology of Taiwan under grant 107-2115-M-006 -007 -MY2.

Department of Mathematics, University of Macau, Macao, China (xqjin@umac.mo). This research was supported in part by the research grant MYRG2016-00077-FST from University of Macau.

(4)

1 Introduction

Let |λ1| ≥ |λ2| ≥ · · · ≥ |λn| ≥ 0 and σ1≥ · · · ≥ σn≥ 0 be the eigenvalues and singular

values of a given n × nrealmatrix A. In [45] Weyl showed that sets of eigenvalues and singular values satisfy the following necessary condition:

k

∏

j=1 |λj| ≤ k

∏

j=1 σj, k= 1, . . . , n − 1, (1.1a) n

∏

j=1 |λj| = n

∏

j=1 σj. (1.1b)

Moreover, Horn [29] proved that condition (1.1), called the Weyl-Horn condition, is also sufficient for constructing triangular matrices with prescribed eigenvalues and singular values. Research interest in inverse eigenvalue and singular value problems can be tracked back to the open problem raised by Higham in [28, Problem 26.3], as follows:

Develop an efficient algorithm for computing a unit upper triangular n × n matrix with the prescribed singular values σ1, . . . , σn, where ∏nj=1σj= 1.

This problem, which was solved by Kosowski and Smoktunowicz [32], leads to the following interesting inverse eigenvalue and singular value problem (IESP):

(IESP) Given two sets of numbers λ = {λ1, . . . , λn} and σ = {σ1, . . . , σn}

satisfying (1.1), find a real n × n matrix with eigenvalues λ and singular values σ .

The following factors make the IESP difficult to solve:

– Often the desired matrices are real. This problem was solved by the authors of [9] with prescribed real eigenvalues and singular values. The method for finding a general real-valued matrix with prescribed complex-conjugate eigenvalues and singular values was also investigated in [33]. In this work, we take an alternative approach to tackle this problem and add further constraints.

– Often the desired matrices are structured. Corresponding to physical applications, the recovered matrices often preserve some common structure such as nonnega-tive entries or predetermined diagonal entries [8, 46]. In this paper, specifically, we offer the condition of the existence of a nonnegative matrix provided that eigenvalues, singular values, and diagonal entries are given. Furthermore, solving the IESP with respect to the diagonal constraint is not enough because entries of the recovered matrices should preserve certain patterns, for example, non-negativity, which correspond to original observations. How to tackle this struc-tured problem is the main thrust of this paper.

The IESP can be regarded as a natural generalization of the inverse eigenvalue problems, which is known for its a wide variety of applications such as the pole as-signment problem [6, 34, 20], applied mechanics [25, 19, 38, 18, 15], and inverse Sturm-Liouville problem [26, 3, 24, 37]. Thus applications of the IESP could be found in wireless communication [39, 17, 43] and quantum information science [21, 30, 46].

(5)

Research results advanced thus far for the IESP do not fully address the above sce-narios. Often, given a set of data, the IESP is studied in parts. That is, there have been extensive investigations of the conditions for the existence of a matrix when the singular values and eigenvalues are provided (i.e., the Weyl-Horn condition [45, 29]), when the singular values and main diagonal entries are provided (i.e., the Sing-Thompson condition [41, 42]), or when the eigenvalues and main diagonal entries are provided (i.e., the Mirsky condition [36]). Also, the above conditions have given rise to numerical approaches, as found in [5, 16, 8, 9, 22, 32, 49].

Our significance in this work is to consider these conditions together. One rela-tively close result is given in [46], where the authors consider a new type of IESP that requires that all three constraints, i.e., eigenvalues, singular values, and diagonal en-tries, be satisfied simultaneously. Theoretically, Wu and Chu generalize the classical Mirsky, Sing-Thompson, and Weyl-Horn conditions and provide one sufficient con-dition for the existence of a matrix with prescribed eigenvalues, singular values, and diagonal entries when n ≥ 3. Numerically, Wu and Chu establish a dynamic system for constructing such a matrix, in which real eigenvalues are given. In this work, we solve an IESP with complex conjugate eigenvalues and with entries fixed at certain locations. Also, we provide the necessary and sufficient condition of the existence of a 2 × 2 nonnegative matrix with prescribed eigenvalues, singular values, and diagonal elements. Note that, in general, the solution of the IESP is not unique or difficult to find once structured requirements are added. To solve an IESP with some specific feature, we combine techniques from differential geometry and for solving nonlinear equations.

We organize this paper as follows. In section 2, we propose the use of the Rie-mannian inexact Newton method for solving an IESP with complex conjugate eigen-values. In section 3, we show that the convergence is quadratic. In section 4, we demonstrate the application of our technique to an IESP with a specific structure that includes nonnegative or predetermined entries to show the robustness and efficiency of our proposed approaches. The concluding remarks and the solvability of the IESP of a 2 × 2 matrix are given in section 5 and the appendix, respectively.

2 Riemannian inexact Newton method

In this section, we explain how the Riemannian inexact Newton method can be ap-plied to the IESP. The problem of optimizing a function on a matrix manifold has received much attention in the scientific and engineering fields due to its peculiarity and capacity. Its applications include, but are not limited to, the study of eigenvalue problems [12, 13, 7, 1, 2, 14, 10, 50, 52, 48, 46, 51], matrix low rank approximation [4, 27], and nonlinear matrix equations [44, 11]. Numerical methods for solving prob-lems involving matrix manifolds rely on interdisciplinary inputs from differential ge-ometry, optimization theory, and gradient flows.

To begin, letO(n) ⊂ Rn×nbe the group of n × n real orthogonal matrices, and let λ = {λ1, . . . , λn} and σ = {σ1, . . . , σn} be the eigenvalues and singular values of an

n× n matrix. We assume without loss of generality that: λ2i−1= αi+βi

√

−1, λ2i= αi−βi

√

(6)

where αi, βi∈ R with βi6= 0 for i = 1, . . . , k, and we define the corresponding block diagonal matrix Λ = diag α1 β1 −β1α1 , . . . , αk βk −βkαk , λ2k+1, . . . , λ2n

and the diagonal matrix

Σ = diag {σ1, . . . , σn} .

Then the IESP is equivalent to finding matrices U, V, Q ∈O(n), and W∈W (n) := {W ∈ Rn×n|Wi, j= 0 if Λi, j6= 0 or i ≥ j, for 1 ≤ i, j ≤ n},

which satisfy the following equation:

F(U,V, Q,W ) = U ΣV>− Q(Λ +W )Q>= 0. (2.1)

Here, we may assume without loss of generality that Q is an identity matrix and simplify Eq. (2.1) as follows:

F(U,V,W ) = U ΣV>− (Λ +W ) = 0. (2.2)

Let X = (U,V,W ) ∈O(n)×O(n)×W (n). Upon using Eq. (2.2), we can see that we might solve the IESP by

finding X ∈O(n) × O(n) × W (n) such that F(X) = 0, (2.3)

where F :O(n)×O(n)×W (n) → Rn×nis continuously differentiable. By making an initial guess, X0, one immediate way to solve Eq. (2.3) is to apply the Newton method

and generate a sequence of iterates by solving

DF(X_k)[∆ Xk] = −F(Xk), (2.4)

for ∆ Xk∈ TXk(O(n) × O(n) × W (n)) and set

Xk+1= RXk(∆ Xk),

where DF(Xk) represents the differential of F at Xkand R is a retraction onO(n) ×

O(n) × W (n). Since Eq. (2.4) is an underdetermined system, it may have more than one solution. Let DF(Xk)∗be the adjoint operator of DF(Xk). In our calculation, we

choose the solution ∆ Xkwith the minimum norm by letting [35, Chap. 6]

∆ Xk= DF(Xk)∗[∆ Zk], (2.5)

where ∆ Zk∈ TF(Xk)(R

n×n_{) is a solution for}

(DF(X_k) ◦ DF(X_k)∗) [∆ Zk] = −F(Xk). (2.6)

Note that the notation ◦ represents the composition of two operators DF(Xk) and

DF(Xk)∗. This implies that the operator DF(Xk) ◦ DF(Xk)∗is symmetric and positive

semidefinite. If, as is the general case, the operator DF(Xk)◦DF(Xk)∗: TF(Xk)(R

n×n_{) →}

(7)

Note that solving for the root of Eq. (2.6) could be unnecessary and computa-tionally time-consuming, and that the linear model given by Eq. (2.6) is large-scale or the resulting iteration Xkis far from the root of condition (2.3) [40]. By analogy

with the classical Newton method [23], we adopt the “inexact” Newton method on Riemannian manifolds, i.e., without solving Eq. (2.6) exactly, we repeatedly apply the conjugate gradient (CG) method to find ∆ Zk∈ TF(Xk)(R

n×n_{), such that:}

k(DF(X_k) ◦ DF(X_k)∗)[∆ Zk] + F(Xk)k ≤ ηkkF(Xk)k, (2.7)

for some constant ηk∈ [0, 1), is satisfied. Then, we update Xkcorresponding to ∆ Zk

until the stopping criterion is satisfied. Here, the notation k · k is the Frobenius norm. Note that in our calculation, the elements in the product space Rn×n× Rn×n_{× R}n×n

are computed using the standard Frobenius inner product:

h(A1, A2, A3), (B1, B2, B3)iF:= hA1, B1i + hA2, B2i + hA3, B3i , (2.8)

where hA, Bi := trace(AB>) for any A, B ∈ Rn×n and the induced norm kX kF =

phX, Xi_F (or, simply, hX , X i and kX k without the risk of confusion) for any X ∈ Rn×n× Rn×n× Rn×n.

Then, the linear mapping DF(X_k) at ∆ Xk= (∆Uk, ∆Vk, ∆Wk) ∈ TXk(O(n)×O(n)×

W (n)) is given by:

DF(Xk)[∆ Xk] = ∆UkΣV_k>+UkΣ ∆V_k>− ∆Wk.

Let DF(Xk)∗: TF(Xk)(R

n×n_{) → T}

Xk(O(n) × O(n) × W (n)) be the adjoint of the

map-ping DF(Xk). The adjoint DF(Xk)∗is determined by the following:

h∆ Zk, DF(Xk)[∆ Xk]i = hDF(Xk)∗[∆ Zk], ∆ Xki

and can be expressed as follows:

DF(Xk)∗[∆ Zk]= (∆Uk, ∆Vk, ∆Wk), where ∆Uk= 1 2(∆ ZkVkΣ > −UkΣVk>∆ Z > kUk), ∆Vk= 1 2(∆ Z > kUkΣ − VkΣ>U_k>∆ ZkVk), ∆Wk= −H ∆ Zk,

with the notation representing the Hadamard product (see [12, 51] for a similar discussion).

There is definitely no guarantee that the application of the inexact Newton method can achieve a sufficient decrease in the size of the nonlinear residual kF(Xk)k. This

provides motivation for deriving an iterate for which the size of the nonlinear residual is decreased. One way to do this is to update the Newton step ∆ Xk obtained from

Eq. (2.5) by choosing θ ∈ [θmin, θmax], with 0 < θmin< θmax< 1, and setting

c

∆ Xk= ∆ Xk, ηˆk=

kF(Xk) + DF(Xk)∆ Xkk

kF(Xk)k

(8)

and ηk= ˆηk. Then, we update ηk← 1 − θ (1 − ηk) and ∆ Xk← 1 − ηk 1 − ˆηk c ∆ Xk, (2.10) while kF(Xk)k − kF(RXk(∆ Xk))k > t(1 − ηk)kF(Xk)k, or, equivalently, kF(RXk(∆ Xk))k < [1 − t(1 − ηk)]kF(Xk)k, (2.11)

for some t ∈ [0, 1) [23]. Let q f (·) denote the mapping that sends a matrix to the Q factor of its QR decomposition with its R factor having strictly positive diagonal ele-ments [1, Example 4.1.3]. Then, for all (ξU, ξV, ξW) ∈ T(U,V,W )(O(n) × O(n) × W (n)),

we can compute the retraction R using the following formula: R(U,V,W )(ξU, ξV, ξW) = (RU(ξU), RV(ξV), RW(ξW)),

where

RU(ξU) = q f (U + ξU), RV(ξV) = q f (V + ξV), RW(ξW) = W + ξW.

We call this the Riemannian inexact Newton backtracking method (RINB) and for-malize this method in Algorithm 1. To choose the parameter θ ∈ [θmin, θmax], we

apply a two-point parabolic model [31, 51] to achieve a sufficient decrease among steps 6 to 9. That is, we use the iteration history to model an approximate minimizer of the following scalar function:

f(λ ) := kF(RXk(λ ∆ Xk))k

2

by defining a parabolic model, as follows:

p(λ ) = f (0) + f0(0)λ + ( f (1) − f (0) − f0(0))λ2,

where f (0) = kF(Xk)k2, f0(0) = 2 hDF(Xk)[∆ Xk], F(Xk)i, and f (1) = kF(RXk(∆ Xk))k

2_.

From (2.7), it can be shown that the function evaluation f0(0) should be negative. Since f0(0) < 0, if p00(λ ) = 2( f (1) − f (0) − f0(0)) > 0, then p(λ ) has its minimum at:

θ = − f

0

(0)

2( f (1) − f (0) − f0₍₀₎₎> 0;

otherwise, if p00(λ ) < 0, we choose θ = θmax. By incorporating two types of selection,

we can choose the following:

θ = min max θmin, − f0₍₀₎ 2( f (1) − f (0) − f0(0)) , θmax .

as the parameter θ in Algorithm 1 [31, 51]. In the next section, we mathematically investigate the convergence analysis of Algorithm 1.

(9)

Algorithm 1: The Riemannian inexact Newton backtracking method [X ] = RINB(σ , X0)

Input: An initial value X0

Output: A numerical solution X satisfying F(X ) = 0 1 begin

2 Let ηmax∈ [0.1), η0= min{ηmax, kF(X0)k}, and t ∈ [0, 1), and 0 < θmin< θmax< 1 be given.

3 repeat

4 Determine ∆ Zkby using the CG method to (2.6) until (2.7) holds.

5 Set ∆ Xk= (DF(Xk))∗∆ Zk, ˆηk=kF(Xk)+DF(XkF(Xk)kk)∆ Xkk, c∆ Xk= ∆ Xk, and ηk= ˆηk.

6 repeat

7 Choose θ ∈ [θmin, θmax].

8 Update ηk← 1 − θ (1 − ηk) and ∆ Xk←1− ˆ1−ηηkk∆ Xck.

9 until (2.11) holds;

10 Set Xk+1= RXk(∆ Xk) and ηk+1= min {ηk, ηmax, kF(Xk+1)k}.

11 Replace k by k + 1. 12 until kF(Xk)k < ε;

13 X= Xk.

14 end

3 Convergence Analysis

By combining the classical inexact Newton method [23] with optimization techniques on matrix manifolds, Algorithm 1 provides a way to solve the IESP. However, we have yet to theoretically discuss the convergence analysis of Algorithm 1. In this sec-tion, we provide a theoretical foundation for the RINB method, and show that this RINB method converges globally and finally converges quadratically when Algo-rithm 1 does not terminate prematurely. We address this phenomenon in the follow-ing:

Lemma 3.1 Algorithm 1 does not break down at some Xkif and only if F(Xk) 6= 0

and the inverse of DF(Xk) ◦ DF(Xk)∗exists.

Next, we provide an upper bound for the approximate solution c∆ Xk in

Algo-rithm 1.

Theorem 3.1 Let ∆ Zk∈ TF(Xk)(R

n×n_{) be a solution that satisfies condition (2.7) and}

c ∆ Xk= DF(Xk)∗[∆ Zk]. Then, (a) k c∆ Xkk ≤ (1 + ˆηk)kDF(Xk)†kkF(Xk)k, (3.1a) (b) kσk(η)k ≤ 1 + ηmax 1 − ηmax (1 − η)kDF(Xk)†kkF(Xk)k, (3.1b)

where ˆηkis defined in Eq.(2.9), and σkis the backtracking curve used in Algorithm 1,

which is defined by the following: σk(η) =

1 − η 1 − ˆηk

c ∆ Xk

(10)

with ˆηk≤ η ≤ 1, and

kDF(Xk)†k := max

k∆ Zk=1kDF(Xk) †

[∆ Z]k represents the norm of the pseudoinverse of DF(X_k).

Proof Let r_k= (DF(X_k) ◦ DF(X_k)∗)[∆ Zk] + F(Xk). We see that

k c∆ Xkk ≤ kDF(Xk)∗◦ [DF(Xk) ◦ DF(Xk)∗]−1kkrk− F(Xk)k ≤ (1 + ˆηk)kDF(Xk)†kkF(Xk)k and kσk(η)k = 1 − η 1 − ˆηk kDF(Xk)†(rk− F(Xk))k ≤ 1 + ˆηk 1 − ˆηk (1 − η)kDF(Xk)†kkF(Xk)k ≤ 1 + ηmax 1 − ηmax (1 − η)kDF(Xk)†kkF(Xk)k. u t In our subsequent discussion, we assume that Algorithm 1 does not break down and there is a unique limit point X∗of {Xk}. Since F is continuously differentiable,

we have the following:

kDF(X)†_{k ≤ 2kDF(X}

∗)†k (3.2)

whenever X ∈ B_δ(X∗) for a sufficiently small constant δ > 0. Here, the notation

B_δ(X∗) represents a neighborhood of X∗consisting of all points X such that kX −

X∗k < δ . By condition (3.1), we can show without any difficulty that whenever Xkis

sufficiently close to X∗,

k c∆ Xkk ≤ (1 + ηmax)kDF(X∗)†kkF(Xk)k, (3.3)

kσk(η)k ≤ Γ (1 − η)kF(Xk)k, ηˆk≤ η ≤ 1,

where Γ is a constant independent of k defined by Γ = 21 + ηmax

1 − ηmax

kDF(X∗)†k.

New, we show that the sequence of {F(Xk)} eventually converges to zero.

Theorem 3.2 Assume that Algorithm 1 does not break down. If {X_k} is the sequence generated in Algorithm 1, then

lim

k→∞F(Xk) = 0.

Proof Observe that

kF(X_k)k = kF(RXk−1(∆ Xk−1))k ≤ (1 − t(1 − ηk−1))kF(Xk−1)k ≤ kF(X0)k k−1

∏

j=0 (1 − t(1 − ηj)) ≤ kF(X0)ke −tk−1∑ j=0 (1−ηj) .

Since t > 0 and lim

k→∞ k−1 ∑ j=0 (1 − ηj) = ∞, we have lim k→∞F(Xk) = 0. ut

(11)

In our iteration, we implement the repeat loop among steps 6 to 9 by selecting a sequence {θj}, with θj∈ [θmin, θmax]. For each loop, correspondingly, we let η

(1) k =

ˆ

ηkand ∆ X(1)= c∆ Xk, and for j = 2, . . . , we let

η_k( j)= 1 − θj−1(1 − ηk( j−1)), ∆ X_k( j)= 1 − η ( j) k 1 − ˆηk c ∆ Xk. (3.4)

By induction, then, we can easily show that:

∆ X_k( j)= Θj−1∆ Xck, 1 − η_k( j)= Θj−1(1 − ˆηk), where Θj−1= j−1

∏

`=1 θ`, j≥ 2. (3.5)

That is, the sequence {∆ X_k( j)}jis a strictly decreasing sequence satisfying lim j→∞∆ X

( j) k =

0, and {η_k( j)}jis a sequence satisfying η ( j)

k ≥ ˆηkfor j ≥ 1, and lim j→∞η

( j)

k = 1. Based

on these observations, next, we show that the repeat loop terminates after a finite number of steps.

Theorem 3.3 Let { c∆ Xk} be the sequence generated from Algorithm 1, i.e.,

k(DF(Xk)[ c∆ Xk] + F(Xk)k ≤ ηkkF(Xk)k.

Then, once j is large enough, the sequence{η_k( j)}jsatisfies the following:

kF(Xk) + DF(Xk)[∆ X ( j) k ]k ≤ η ( j) k kF(Xk)k, kF(RXk(∆ X ( j) k ))k ≤ (1 − t(1 − η ( j) k ))kF(Xk)k. (3.6)

Proof Let ˆηkbe defined in Eq. (2.9) with ∆ Xk= c∆ Xk, and εk=(1−t)(1− ˆ_{k c}ηk)kF(Xk)k ∆ Xkk

. Since F is continuously differentiable, for εk> 0, there exists a sufficiently small

δ > 0 such that k∆ X k < δ implies that:

kF(RXk(∆ X )) − F(RXk(0Xk)) − DF(RXk(0Xk))[∆ X ]k ≤ εkk∆ Xk,

where 0Xk is the origin of TXk(O(n) × O(n) × W (n)).

For δ > 0, we let ηmin= max ( ˆ ηk, 1 − (1 − ˆηk)δ k c∆ Xkk ) .

Note that once j is sufficiently large,

η_k( j)− ηmin≥ δ k c∆ Xkk −Θj−1 ! (1 − ˆηk)≥0. (3.7)

(12)

For sufficiently large j, we consider the sequence {∆ X_k( j)}jin Eq. (3.4) with η_k( j)∈

[ηmin, 1). We can see that:

k∆ X_k( j)k = k1 − η ( j) k 1 − ˆηk c ∆ Xkk ≤ 1 − ηmin 1 − ˆηk k c∆ Xkk ≤ δ .

This implies that: kF(Xk) + DF(Xk)[∆ X ( j) k ]k ≤ F(X_k) + DF(X_k) 1 − η ( j) k 1 − ˆηk c ∆ Xk ! ≤ η_k( j)− ˆηk 1 − ˆηk F(Xk) + 1 − η_k( j) 1 − ˆηk DF(Xk)[ c∆ Xk] + F(Xk) ≤ η ( j) k − ˆηk 1 − ˆηk kF(Xk)k + 1 − η_k( j) 1 − ˆηk ˆ ηkkF(Xk)k = η_k( j)kF(X_k)k, and kF(RXk(∆ X ( j) k ))k= kF(RXk(∆ X ( j) k ) − F(RXk(0Xk)) − DF(RXk(0Xk))[∆ X ( j) k ]k +kF(Xk) + DF(Xk)[∆ X ( j) k ]k ≤εkk∆ X ( j) k k + η ( j) k kF(Xk)k = (1 − t)(1 − ˆηk)kF(Xk)k k c∆ Xkk 1 − η_k( j) 1 − ˆηk c ∆ Xk + η_k( j)kF(Xk)k = (1 − t(1 − η_k( j)))kF(Xk)k. u t From the proof of Theorem 3.3, we can see that for each k, the repeat loop for the backtracking line search will terminate in a finite number of steps once condition (3.7) is satisfied. Moreover, Theorem 3.2 and condition (3.3) imply the following:

lim

k→∞k c∆ Xkk = 0.

That is, if k is sufficient large, i.e., k c∆ Xkk is small enough, then from the proof of

Theorem 3.3 we see that condition (2.11) is always satisfied, i.e., ηk= ˆηk for all

sufficient large k.

To show that Algorithm 1 is a globally convergent algorithm, we have one ad-ditional requirement for the retraction RX, i.e., there exist ν > 0 and δν > 0 such

that:

ν k∆ X k ≥ dist(RX(∆ X ), X ), (3.8)

for all X ∈O(n) × O(n) × W (n) and for all ∆X ∈ TX(O(n) × O(n) × W (n)) with

k∆ Xk ≤ δν[1]. Here “dist(·, ·)” represents the Riemannian distance onO(n)×O(n)×

W (n). Under this assumption, our next theorem shows the global convergence prop-erty of Algorithm 1. We have borrowed the strategy for this proof from that used in [23, Theorem 3.5] to prove the nonlinear matrix equation.

(13)

Theorem 3.4 Assume that Algorithm 1 does not break down. Let X∗be a limit point

of {Xk}. Then Xk converges to X∗ and F(X∗) = 0. Moreover, Xk converges to X∗

quadratically whenever Xkis sufficiently close to X∗.

Proof Suppose Xk does not converge to X∗. This implies that there exist two

se-quences of numbers {kj} and {`j} for which:

Xkj ∈ Nδ / j(X∗),

X_k_j+`j 6∈ Nδ(X∗),

Xkj+i∈ Nδ(X∗), if i = 1, . . . , `j−1

kj+ `j≤ kj+1.

From Theorem 3.3, we see that the repeat loop among steps 6 to 9 of Algorithm 1 terminates in finite steps. For each k, let mkbe the smallest number such that

condi-tion (3.6) is satisfied, i.e., ∆ Xk= Θmk∆ Xckand ηk= 1 −Θmk(1 − ˆηk) with Θmk being

defined in Eq. (3.5). It follows from condition (3.1b) that: k∆ Xkk ≤ 2Θmk

1 + ηmax

1 − ηmax

(1 − ηk)kDF(X∗)†kkF(Xk)k, (3.9)

for a sufficiently small δ and Xk∈ Bδ(X∗), so that condition (3.2) is satisfied. Let

Γmk= 2Θmk

1 + ηmax

1 − ηmax

kDF(X∗)†k.

According to condition (3.8), there exist ν > 0 and δν> 0 such that:

ν k∆ X k ≥ dist (RX(∆ X ), X ) ,

when k∆ X k ≤ δν. Since F(Xk) approaches zero as k approaches infinity, for δν,

con-dition (3.9) implies that there exists a sufficiently large k such that: ν k∆ Xkk ≥ dist RXk(∆ Xk), Xk

(3.10) is satisfied whenever k∆ Xkk ≤ δν.

Then for a sufficiently large j, we can see from conditions (3.9) and (3.10) that: δ 2 ≤ dist(Xkj+`j, Xkj) ≤ kj+`j−1

∑

k=kj dist(Xk+1, Xk) = kj+`j−1

∑

k=kj dist(RXk(∆ Xk), Xk) ≤ kj+`j−1

∑

k=kj ν k∆ Xkk ≤ kj+`j−1

∑

k=kj νΓmk(1 − ηk)kF(Xk)k ≤ kj+`j−1

∑

k=kj νΓmk t (kF(Xk)k − kF(Xk+1)k) ≤ νΓmk t kF(Xkj)k − kF(Xkj+1)k .

(14)

This is a contraction, since Theorem 3.2 implies that F(Xkj) converges to zero as j

approaches infinity and Γmk is bounded. Thus, Xkconverges to X∗, and immediately,

we have F(X∗) = 0. This completes the proof of the first part.

To show that Xkconverges to X∗quadratically once Xkis sufficiently close to X∗,

we let C1and C2be two numbers satisfying the following:

kF(Xk+1) − F(Xk) − DF(Xk)[∆ Xk]k ≤ C1k∆ Xkk2,

kF(Xk)k ≤ C2dist(Xk, X∗),

for a sufficiently large k. The above assumptions are true since F is second differen-tiable and F(X∗) = 0. We can also observe that:

kF(Xk+1)k ≤ kF(Xk+1) − F(Xk) − DF(Xk)[∆ Xk]k + kF(Xk) + DF(Xk)[∆ Xk]k ≤ C1k∆ Xkk2+ ˆηkkF(Xk)k ≤ C1(ΓmkkF(Xk)k) 2_{+ kF(X} k)k2 ≤ C1Γ2C₂2+C₂2 dist(Xk, X∗)2, (3.11) where Γ = 2 1 + ηmax 1 − ηmax kDF(X∗)†k.

Since Xk converges to X∗as k converges to infinity, for a sufficiently large k, it

follows from conditions (3.9), (3.10), (3.6), and (3.11) that: dist(Xk+1, X∗) = lim p→∞dist(Xk+1, Xp) ≤ ∞

∑

s=k dist Xs+1, RXs+1(∆ Xs+1) ≤ ∞

∑

s=k ν k∆ Xs+1k ≤ ∞

∑

s=k νΓms+1(1 + ηmax)kF(Xs+1)k ≤ νΓ (1 + ηmax) ∞

∑

j=0 (1 − t(1 − ηmax))jkF(Xk+1)k ≤ Cdist(Xk, X∗)2,

for some constant C = νΓ (1 + ηmax) C1Γ

2_C2 2+C22

t(1 − ηmax)

. ut

It is true that we might assume without loss of generality that the inverse of DF(Xk) ◦ DF(Xk)∗always exists numerically. However, once DF(Xk) ◦ DF(Xk)∗is

ill-conditioned or (nearly) singular, we choose an operator Ek= σkidTF_(Xk), where σk

is a constant and idTF_(Xk) is an identity operator on TF(Xk)(R

n×n_{) to make DF(X}

k) ◦

DF(X_k)∗+σkidTF_(Xk)well-conditioned or nonsingular. In the calculation, this replaces

the calculation in Eq. (2.6) with the following equation:

(DF(X_k) ◦ DF(X_k)∗+ σkidTF_(Xk))[∆ Zk] = −F(Xk).

That is, Algorithm 1 can be modified to fit in this case by replacing the satisfaction of condition (2.7) with the following two conditions:

k(DF(Xk) ◦ DF(Xk)∗+ σkidTF_(Xk))[∆ Zk]k ≤ ηkkF(Xk)k, (3.12a)

(15)

where σk:= min {σmax, kF(Xk)k} is a selected perturbation determined by the

param-eter σmaxand kF(Xk)k. Of course, we can provide the proof of the quadratic

conver-gence under condition (3.12) without any difficulty (see [51] for a similar discussion). Thus, we ignore the proof here. However, we note that even if a selected perturba-tion is applied to an ill-condiperturba-tioned problem, the linear operator DF(Xk) ◦ DF(Xk)∗+

σkidTF_(Xk) in condition (3.12a) might become nearly singular or ill-conditioned once

σkis small enough. This will prevent the iteration in the CG method from converging

in fewer than n2steps, and cause the value of f0(0) to not be negative. This possibility suggests that we apply Algorithm 1 without performing any perturbation in our nu-merical experiments. If the CG method cannot terminate within n2_{iterations, it may}

be necessary to compute a new approximated solution ∆ Zkby selecting a new initial

value for X0.

4 Numerical Experiments

Note that the iteration of Algorithm 1 will be trapped without convergence to a solu-tion if the IESP is unsolvable. As such, in our numerical experiments, we assume the existence of a solution of an IESP solution beforehand by generating sets of eigen-values and singular eigen-values from a series of randomly generated matrices. For a 2 × 2 case, it is certain that Theorem A.3 in the appendix provides an alternative way to generate testing matrices. However, for general n × n matrices, the condition of the solvability of the IESP with some particular structure remains unknown and mer-its further investigation. In this section, we show how Algorithm 1 can be applied to solve an IESP with a particular structure. We note that we performed all of the computations in this work in MATLAB version 2016a on a desktop with a 4.2 GHZ Intel Core i7 processor and 32 GB of main memory. For our tests, we set ηmax= 0.9,

θmin= 0.1, θmax= 0.9, t = 10−4, and ε = 10−10. Also, in our computation, we

em-phasize two things. First, once the CG method computed in Algorithm 1 cannot be terminated within n2iterations, restart Algorithm 1 with a different initial value X0.

Second, due to the rounding errors in numerical computation, care must be taken in the selection of ηkso that the upper bound ηkkF(Xk)k in condition (2.7) is not too

small to cause the CG method abnormal. To this end, in our experiments, we use the condition

max{ηkkF(Xk)k, 10−12},

instead of ηkkF(Xk)k. The implementations of the Algorithm 1 are available online,

say, http://myweb.ncku.edu.tw/~mhlin/Bitcodes.zip.

Example 4.1 To demonstrate the capacity of our approach for solving problems that are relatively large, we randomly generate a set of eigenvalues and a set of singular values of different size, say, n = 20, 60, 100, 150, 200, 500, and 700 from matrices given by the MATLAB command:

(16)

For each size, we perform 10 experiments. To illustrate the elasticity of our approach, we randomly generate the initial value X0= (U0,V0,W0) in the following way:

W0= triu(randn(n)), W0(find(Λ)) = 0, and [U0,tmp,V0] = svd(Λ +W0).

In Table 4.1, we show the average residual value (Residual), the average final error (Error), as defined by:

final error = kλ (Anew) − λ k2+ kσ (Anew) − σ k2,

the average number of iterations within the CG method (CGIt)], the average number of iterations within the inexact Newton method (INMIt)], and the average elapsed time (Time), as performed by our algorithm. In Table 4.1, we can see that the elapsed time and the average number of iterations within the CG method increase dramat-ically as the size of the matrices increases. This can be explained by the fact that the number of degrees of freedom of the problem increases significantly. Thus, the number of the iterations required by the CG method and the required computed time increase correspondingly. However, it is interesting to see that the required number of iterations within the inexact Newton method remains almost the same for matrices of different sizes. One way to speed up the entire process of iterations is to transform the problem (2.6) into a form that is more suitable for the CG method, for example, apply the CG method with a preselected preconditioner. Still, this selection of the preconditioner requires further investigation.

Table 4.1 Comparison of the required CGIt], INMIt], Residual, Error values, and Time for solving the IESP by Algorithm 1.

n CGIt] INMIt] Residual Error Time 20 208 9.4 5.54 × 10−12 9.65 × 10−13 2.47 × 10−2 60 740 10 8.13 × 10−12 7.23 × 10−13 4.11 × 10−1 100 1231 10.4 1.06 × 10−12 9.74 × 10−14 2.22 150 1773 10.1 1.01 × 10−12 1.06 × 10−13 6.82 200 1939 10.5 1.20 × 10−12 1.49 × 10−13 19.3 500 6070 10.6 1.47 × 10−12 4.12 × 10−13 665 700 8905 10.6 5.42 × 10−12 7.24 × 10−13 2465

Example 4.2 In this example, we use Algorithm 1 to construct a nonnegative matrix with prescribed eigenvalues and singular values and a specific structure. We specify this IESP and call it the IESP with desired entries (DIESP). The DIESP can be defined as follows.

(DIESP) Given a subsetI = {(it, jt)}`t=1 with double subscripts, a set of

real numbers K = {kt}`t=1, a set of n complex numbers {λi}ni=1,

satisfy-ing {λi}ni=1= { ¯λi}ni=1, and a set of n nonnegative numbers {σi}ni=1, find a

nonnegative n × n matrix A that has eigenvalues λ1, . . . , λn, singular values

(17)

Note that once it= jt= t for t = 1, . . . , n, we investigate a numerical approach for

solving the IESP with prescribed diagonal entries. As far as we know, the research result close to this problem is only available in [46]. However, for a general structure, no research has been conducted to implement this investigation. To solve the DIESP, our first step is to obtain a real matrix A with prescribed eigenvalues and singular values. Our second step is to derive entries of Q>AQ, where Q ∈O(n), that satisfy the nonnegative property and desired values determined by the setsI and K . We solve the first step in the same manner as in Example 4.1, but for the second step, we consider the following two setsL1andL2, which are defined by:

L1= {A ∈ Rm×n| Ait, jt = kt, for 1 ≤ t ≤ `; otherwise Ai, j= 0},

L2= {A ∈ Rm×n| Ai, j= 0, for 1 ≤ i, j ≤ n and (i, j) ∈I },

and then solve the following problem:

find P ∈L2and Q ∈O(n) such that H(P,Q) = ˆA + P P − QAQ>= 0, (4.1)

with ˆA∈L1. Let [A, B] := AB − BA denote the Lie bracket notation. It follows from

direct computation that the corresponding differential DH and its adjoint DH∗have the following form [51]:

DH(P, Q)[(∆ P, ∆ Q)] = 2P ∆ P + [QAQ>, ∆ QQ>], DH(P, Q)∗[∆ Z] = 2P ∆ Z,1 2([QAQ > , ∆ Z>] + [QA>Q>, ∆ Z])Q , and, for all (ξP, ξQ) ∈ T(P,Q)(L2×O(n)), we can compute the retraction R using the

following formula:

R(P, Q) = (RP(ξP), RQ(ξQ)),

where

RP(ξP) = P + ξP, RQ(ξQ) = q f (Q + ξQ).

For these experiments, we randomly generate nonnegative matrices 20 × 20 in size by the MATLAB command “A = rand(20)” to provide the desired eigenvalues, singular values, and diagonal entries, i.e., to solve the DIESP with the specified diag-onal entries. We record the final error, as given by the following formula:

final error = kλ (Anew) − λ k2+ kσ (Anew) − σ k2+ k(Anew)it, jt− ktk2.

After randomly choosing 10 different matrices, Table 4.2 shows our results with the intervals (Interval) containing all of the residual values and final errors, and their corresponding average values (Average). These results provide sufficient evidence that Algorithm 1 can be applied to solve the DIESP with high accuracy.

Although Example 4.2 considers examples with a nonnegative structure, we em-phasize that Algorithm 1 can work with entries that are not limited to being non-negative. That is, to solve the IESP without nonnegative constraints but with another specific structure, Algorithm 1 can fit perfectly well by replacing H(P, Q) in prob-lem (4.1) with

G(S, Q) := ˆA+ S − QAQ>, where ˆA∈L1, S ∈L2and Q ∈O(n).

(18)

Table 4.2 Records of final errors and residual values for solving the DIESP by Algorithm 1. Interval Average final errors [7.27 × 10−13, 1.21 × 10−11] 2.91 × 10−12 residual values [7.77 × 10−13_{, 4.93 × 10}−12_] _{1.85 × 10}−12

5 Conclusions

In this paper, we apply the Riemannian inexact Newton method to solve an initially complicated and challenging IESP. We provide a thorough analysis of the entire iter-ative processes and show that this algorithm converges globally and quadratically to the desired solution. We must emphasize that our theoretical discussion and numeri-cal implementations can also be extended to solve an IESP with a particular structure such as desired diagonal entries and a matrix whose entries are nonnegative. This capacity can be observed in our numerical experiments. It should be emphasized that this research is the first to provide a unified and effective means to solve the IESP with or without a particular structure.

However, the numerical stability for extremely ill-conditioned problems is a case that we should pay attention to, though reselecting the initial values could be a strat-egy to get rid of this difficulty. Another way to tackle this difficulty is to select a good preconditioner. But, the operator encountered in our algorithm is nonlinear and high-dimensional. Thus, the selection of the preconditioner could involve the study of tensor analysis, where further research is needed.

Theoretically determining the sufficient and necessary condition for solving IESPs of any specific structure, including a stochastic, Toeplitz, or Hankel structure, is chal-lenging and interesting. In the appendix, we provide the solvability condition of the IESP with real or nonnegative matrices of size 2 × 2 real/nonnegative matrices, while the desired eigenvalues, singular values, and main diagonal entries are given. We hope that this discussion can motivate a further discussion shortly.

A Appendix

A.1 The solvability of the IESP of a 2 × 2 matrix

For the IESP, the authors in [46] use a geometric argument to investigate a necessary and sufficient condi-tion for the existence of a 2 × 2 real matrix with prescribed diagonal entries. This argument also leads to a sufficient algebraic but not necessary condition for the construction of a 2 × 2 real matrix. In this appendix, the algebraic condition under which a 2 × 2 real matrix or even nonnegative matrix can be constructed in closed form, given its eigenvalue, singular values, and main diagonal entries. To do so, we must have the following results. The first result, the so-called Mirsky condition, provides the classical relationship between the eigenvalues λ = {λ1, . . . , λn} and the diagonal entries d = {d1, . . . , dn}.

Theorem A.1 [[36], Mirsky condition]. There exists a real matrix A ∈ Rn×n having eigenvalues λ = {λ1, . . . , λn} and main diagonal entries d = {d1, . . . , dn}, that are possibly in different order, if and only if

n

∑

i=1 λi= n

∑

i=1 di. (A.1)

(19)

The second result provides the relationship between the singular values σ and main diagonal entries dof a 2 × 2 nonnegative matrix.

Theorem A.2 [[47], Theorem 2.1]. There exists a nonnegative matrix A = d1 b c d2

∈ R2×2_{having the}

singular values σ1≥ σ2and main diagonal entries d1≥ d2, with renumbering if necessary, if and only if

σ1+ σ2≥ d1+ d2, σ1− σ2≥ d1− d2, if bc− d1d2≤ 0, (A.2a)

σ1− σ2≥ d1+ d2, if bc− d1d2> 0. (A.2b)

In particular, entries from matrix A can be relaxed to real numbers, and condition (A.2) is also true for the construction of a 2 × 2 real matrix. The proof is almost identical to that in [47, Lemma 2.1]. The major change is the substitution of nonnegative entries for real entries. Thus, we skip its proof here.

Theorem A.3 There exists a real matrix A = d1 b c d2

∈ R2×2_{having singular values σ}

1≥ σ2and main

diagonal entries d1≥ d2, with renumbering if necessary, if and only if

σ1+ σ2≥ d1+ d2, σ1− σ2≥ d1− d2, if bc− d1d2≤ 0,

σ1− σ2≥ d1+ d2, if bc− d1d2> 0.

Now we have the condition of the existence of a 2 × 2 matrix provided with eigenvalues and main diagonal entries, or singular values and main diagonal entries. The next theorem, unsolved in [47], deals with the case in which the three constraints—eigenvalues, singular values, and main diagonal entries—are of simultaneous concern.

Theorem A.4 There exists a real matrix A = d1 b c d2

∈ R2×2_{having eigenvalues}_|λ

1| ≥ |λ2|, singular

values σ1≥ σ2, and main diagonal entries d1≥ d2, with renumbering if necessary, if and only if

λ1+ λ2= d1+ d2, σ1≥ |λ1|, |λ1λ2| = σ1σ2, (A.4)

and

σ1+ σ2≥ d1+ d2, σ1− σ2≥ d1− d2, if bc− d1d2≤ 0, (A.5a)

σ1− σ2≥ d1+ d2, if bc− d1d2> 0. (A.5b)

Proof Assume that conditions (A.4) and (A.5) are satisfied. Following from the Weyl-Horn and Mirsky conditions, we know that for any 2 × 2 matrix, its eigenvalues, singular values, and diagonal entries must satisfy condition (A.4). Thus, Theorem A.3 implies that once condition (A.5) is satisfied, it suffices to say that there exists a 2 × 2 real matrix.

On the other hand, the sufficient condition follows directly from the Weyl-Horn condition (1.1), the Mirsky condition (A.1), and Theorem A.3. This completes the proof. ut Since the solvability conditions of Theorem A.2 and Theorem A.3 are equivalent, we can see that the solvability condition in Theorem A.4 can be confined to be the necessary and sufficient condition for the existence of a nonnegative 2 × 2 matrix. We summarize this result as follows.

Corollary A.1 There exists a nonnegative matrix A = d1 b c d2

∈ R2×2_{having eigenvalues}_|λ 1| ≥ |λ2|,

singular values σ1≥ σ2, and main diagonal entries d1≥ d2, with renumbering if necessary, if and only if

λ1+ λ2= d1+ d2, σ1≥ |λ1|, |λ1λ2| = σ1σ2,

and

σ1+ σ2≥ d1+ d2, σ1− σ2≥ d1− d2, if bc− d1d2≤ 0,

σ1− σ2≥ d1+ d2, if bc− d1d2> 0.

Note that conditions (A.4) and (A.5) cannot be directly generalized to higher dimensional cases. The authors in [47] present the necessary and sufficient condition of the existence of a real matrix with a size greater than 2 and having prescribed eigenvalues, singular values, and main diagonal entries. However, given eigenvalues, singular values, and main diagonal entries, no study has yet demonstrated the construc-tion of a nonnegative matrix with a size greater than 2 × 2. This difficulty can be tackled by the use of our numerical computations.

(20)

Acknowledgment

The authors wish to thank Prof. Michiel E. Hochstenbach for his highly valuable comments. They also thank Prof. Zheng-Jian Bai and Dr. Zhi Zhao for helpful discussions.

References

1. P.-A. Absil, R. Mahony, and R. Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ, 2008.

2. C. G. Baker, P.-A. Absil, and K. A. Gallivan. An implicit Riemannian trust-region method for the symmetric generalized eigenproblem. In Computational Science – ICCS 2006, LNCS. Springer, 2006. 3. D. L. Boley and G. H. Golub. A survey of matrix inverse eigenvalue problems. Inverse Problems,

3(4):595–622, 1987.

4. S. Boyd and L. Xiao. Least-squares covariance matrix adjustment. SIAM J. Matrix Anal. Appl., 27(2):532–546, 2005.

5. N. N. Chan and K. H. Li. Diagonal elements and eigenvalues of a real symmetric matrix. J. Math. Anal. Appl., 91(2):562–566, 1983.

6. E. K. Chu and B. N. Datta. Numerically robust pole assignment for second-order systems. Internat. J. Control, 64(6):1113–1127, 1996.

7. M. T. Chu. Numerical methods for inverse singular value problems. SIAM J. Numer. Anal., 29(3):885– 903, 1992.

8. M. T. Chu. On constructing matrices with prescribed singular values and diagonal elements. Linear Algebra Appl., 288(1-3):11–22, 1999.

9. M. T. Chu. A fast recursive algorithm for constructing matrices with prescribed eigenvalues and singular values. SIAM J. Numer. Anal., 37(3):1004–1020, 2000.

10. M. T. Chu. Linear algebra algorithms as dynamical systems. Acta Numer., 17:1–86, 2008.

11. M. T. Chu. On the first degree Fej´er-Riesz factorization and its applications to X + A∗X−1A= Q. Linear Algebra Appl., 489:123–143, 2016.

12. M. T. Chu and K. R. Driessel. The projected gradient method for least squares matrix approximations with spectral constraints. SIAM J. Numer. Anal., 27(4):1050–1060, 1990.

13. M. T. Chu and K. R. Driessel. Constructing symmetric nonnegative matrices with prescribed eigen-values by differential equations. SIAM J. Math. Anal., 22(5):1372–1387, 1991.

14. M. T. Chu and G. H. Golub. Inverse Eigenvalue Problems: Theory, Algorithms, and Applications. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 2005. 15. M. T. Chu, W.-W. Lin, and S.-F. Xu. Updating quadratic models with no spillover effect on

unmea-sured spectral data. Inverse Problems, 23(1):243–256, 2007.

16. M. T. Chu and J. W. Wright. The educational testing problem revisited. IMA Journal of Numerical Analysis, 15(1):141–160, 1995.

17. P. Cotae and M. Aguirre. On the construction of the unit tight frames in code division multiple access systems under total squared correlation criterion. AEU - International Journal of Electronics and Communications, 60(10):724 – 734, 2006.

18. B. N. Datta. Finite-element model updating, eigenstructure assignment and eigenvalue embedding techniques for vibrating systems. Mech. Sys. Signal Processing, 16(1):83 – 96, 2002.

19. B. N. Datta, S. Elhay, Y. M. Ram, and D. R. Sarkissian. Partial eigenstructure assignment for the quadratic pencil. J. Sound Vibration, 230(1):101–110, 2000.

20. B. N. Datta, W.-W. Lin, and J.-N. Wang. Robust partial pole assignment for vibrating systems with aerodynamic effects. IEEE Trans. Automat. Control, 51(12):1979–1984, 2006.

21. A. S. Deakin and T. M. Luke. On the inverse eigenvalue problem for matrices (atomic corrections). Journal of Physics A: Mathematical and General, 25(3):635, 1992.

22. I. S. Dhillon, R. W. Heath, Jr., M. A. Sustik, and J. A. Tropp. Generalized finite algorithms for constructing Hermitian matrices with prescribed diagonal and spectrum. SIAM J. Matrix Anal. Appl., 27(1):61–71, 2005.

23. S. C. Eisenstat and H. F. Walker. Globally convergent inexact Newton methods. SIAM J. Optim., 4(2):393–422, 1994.

24. G. M. L. Gladwell. Inverse Problems in Vibration, volume 119 of Solid Mechanics and its Applica-tions. Kluwer Academic Publishers, Dordrecht, second edition, 2004.

(21)

25. I. Gohberg, P. Lancaster, and L. Rodman. Matrix Polynomials. Society for Industrial and Applied Mathematics, 2009.

26. G. H. Golub. Some modified matrix eigenvalue problems. SIAM Review, 15(2):318–334, 1973. 27. I. Grubiˇsi´c and R. Pietersz. Efficient rank reduction of correlation matrices. Linear Algebra Appl.,

422(2-3):629–653, 2007.

28. N. J. Higham. Accuracy and Stability of Numerical Algorithms. SIAM, Philadelphia, PA, 1996. 29. A. Horn. On the eigenvalues of a matrix with prescribed singular values. Proc. Amer. Math. Soc.,

5:4–7, 1954.

30. K. Jacek. Inverse problems in quantum chemistry. International Journal of Quantum Chemistry, 109(11):2456–2463, 2009.

31. C. T. Kelley. Iterative Methods for Linear and Nonlinear equations, volume 16 of Frontiers in Applied Mathematics. SIAM, Philadelphia, PA, 1995.

32. P. Kosowski and A. Smoktunowicz. On constructing unit triangular matrices with prescribed singular values. Computing, 64(3):279–285, 2000.

33. C.-K. Li and R. Mathias. Construction of matrices with prescribed singular values and eigenvalues. BIT Numerical Mathematics, 41(1):115–126, 2001.

34. N. Li. A matrix inverse eigenvalue problem and its application. Linear Algebra Appl., 266:143–152, 1997.

35. D. G. Luenberger. Optimization by Vector Space Methods. John Wiley & Sons, Inc., New York-London-Sydney, 1969.

36. L. Mirsky. Matrices with prescribed characteristic roots and diagonal elements. J. London Math. Soc., 33:14–21, 1958.

37. M. M¨oller and V. Pivovarchik. Inverse Sturm–Liouville Problems. Springer International Publishing, Cham, 2015.

38. N. K. Nichols and J. Kautsky. Robust eigenstructure assignment in quadratic matrix polynomials: nonsingular case. SIAM J. Matrix Anal. Appl., 23(1):77–102, 2001.

39. R. Rao and S. Dianat. Basics of Code Division Multiple Access (CDMA). SPIE tutorial texts. SPIE Press, 2005.

40. J. P. Simonis. Inexact Newton methods applied to under-determined systems. Ph.D. dissertation, Worcester Polytechnic Institute, Worcester, MA, 2006.

41. F. Y. Sing. Some results on matrices with prescribed diagonal elements and singular values. Canad. Math. Bull., 19(1):89–92, 1976.

42. R. C. Thompson. Singular values, diagonal elements, and convexity. SIAM J. Appl. Math., 32(1):39– 63, 1977.

43. J. A. Tropp, I. S. Dhillon, and R. W. Heath, Jr. Finite-step algorithms for constructing optimal CDMA signature sequences. IEEE Trans. Inform. Theory, 50(11):2916–2921, 2004.

44. L. Wang, M. T. Chu, and Y. Bo. A computational framework of gradient flows for general linear matrix equations. Numer. Algorithms, 68(1):121–141, 2015.

45. H. Weyl. Inequalities between the two kinds of eigenvalues of a linear transformation. Proc. Nat. Acad. Sci. U. S. A., 35:408–411, 1949.

46. S.-J. Wu and M. T. Chu. Solving an inverse eigenvalue problem with triple constraints on eigenvalues, singular values, and diagonal elements. Inverse Problems, 33(8):085003, 21, 2017.

47. S.-J. Wu and M. M. Lin. Numerical methods for solving nonnegative inverse singular value problems with prescribed structure. Inverse Problems, 30(5):055008, 14, 2014.

48. T.-T. Yao, Z.-J. Bai, Z. Zhao, and W.-K. Ching. A Riemannian Fletcher-Reeves conjugate gradient method for doubly stochastic inverse eigenvalue problems. SIAM J. Matrix Anal. Appl., 37(1):215– 234, 2016.

49. H. Zha and Z. Zhang. A note on constructing a symmetric matrix with specified diagonal entries and eigenvalues. BIT, 35(3):448–451, 1995.

50. Z. Zhao, Z.-J. Bai, and X.-Q. Jin. A Riemannian Newton algorithm for nonlinear eigenvalue problems. SIAM J. Matrix Anal. Appl., 36(2):752–774, 2015.

51. Z. Zhao, Z.-J. Bai, and X.-Q. Jin. A Riemannian inexact Newton-CG method for constructing a nonnegative matrix with prescribed realizable spectrum, Numer. Math., published online, 2018. https://doi.org/10.1007/s00211-018-0982-2.

52. Z. Zhao, X.-Q. Jin, and Z.-J. Bai. A geometric nonlinear conjugate gradient method for stochastic inverse eigenvalue problems. SIAM J. Numer. Anal., 54(4):2015–2035, 2016.

(22)

1

科技部補助專題研究計畫執行國際合作與移地研究心得報告

本計畫本人原預計與澳門大學(University of Macau) Xiao-Qing Jin和Haiwei Sun教授進行學術討

論與合作，但由於在近期被通知獲得科技部哥倫布計畫的補助且得提前中斷該執行計畫，故尚未能使

用移地研究經費。還請審查人員見諒。

(23)

1

科技部補助專題研究計畫出席國際學術會議心得報告

本計畫本人原擬出席 SIAM 所舉辦的相關活動，但因為研究的需求改訂參與今年七月在西班牙

所舉辦的 ICIAM 學術活動。但這活動的舉辦日期為 2019/7/15-2019/7/19，由於在近期被通知獲得

科技部哥倫布計畫的補助並得提前中斷該執行計畫，故尚未使用出席國際會議經費。還請審查人員

見諒。

(24)

107年度專題研究計畫成果彙整表

計畫主持人：

林敏雄

計畫編號：

107-2115-M-006-007-MY2

計畫名稱：