Generation of test matrices with specified eigenvalues using floating-point arithmetic

This paper concerns test matrices for numerical linear algebra using an error-free transformation of floating-point arithmetic. For specified eigenvalues given by a user, we propose methods of generating a matrix whose eigenvalues are exactly known based on, for example, Schur or Jordan normal form and a block diagonal form. It is also possible to produce a real matrix with specified complex eigenvalues. Such test matrices with exactly known eigenvalues are useful for numerical algorithms in checking the accuracy of computed results. In particular, exact errors of eigenvalues can be monitored. To generate test matrices, we first propose an error-free transformation for the product of three matrices YSX. We approximate S by S′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${S^{\prime }}$\end{document} to compute YS′X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}${YS^{\prime }X}$\end{document} without a rounding error. Next, the error-free transformation is applied to the generation of test matrices with exactly known eigenvalues. Note that the exactly known eigenvalues of the constructed matrix may differ from the anticipated given eigenvalues. Finally, numerical examples are introduced in checking the accuracy of numerical computations for symmetric and unsymmetric eigenvalue problems.

methods for problems in numerical linear algebra. There are several matrices whose exact eigenvalues are known, for example, Circular, Clement, Frank, and tridiagonal Toeplitz matrices. Test matrices are well summarized in [4,Chapter 28]. This paper aims to generate a matrix with exactly known eigenvalues based on some standard forms, for example, the Schur or the Jordan normal forms, or a block diagonal form. Let X ∈ R n×n be a non-singular matrix. For a matrix S ∈ R n×n , we consider a form: If (1) is the Schur normal form, X is an orthogonal matrix, and S is an upper triangular matrix. For the Jordan normal form, S ∈ R n×n is a matrix whose diagonal blocks are Jordan matrices. It indicates that S is an upper bidiagonal matrix and the sub-diagonal elements s i,i+1 are either 1 or 0. For both normal forms, the diagonal elements s ii (1 ≤ i ≤ n) are eigenvalues of A. Let F be a set of floating-point numbers as defined in IEEE 754 [5]. If we use floating-point numbers and floating-point arithmetic, then there are following problems to get a matrix with exactly known eigenvalues using (1).
Problem 1 Even if all elements in a non-singular matrix X can be represented by floating-point numbers (X ∈ F n×n ), elements in X −1 may not be represented by floating-point numbers (X −1 ∈ F n×n ). Problem 2 Even if X, S, X −1 ∈ F n×n , X −1 SX ∈ F n×n is not always satisfied.
We first prepare a diagonal matrix S and a non-singular matrix X, next obtain the inverse matrix X −1 using exact computations such as symbolic or rational arithmetic. Finally, we compute X −1 SX using exact computations, and the result X −1 SX is rounded to a floating-point matrix A. One may think s ii (1 ≤ i ≤ n) are nearly eigenvalues of A, and these are useful as numerical tests. However, the Bauer-Fike theorem [2] states min λ |λ − μ| ≤ cond(X) E , (2) where λ and μ are eigenvalues of A and A + E = X −1 SX, respectively. Therefore, s ii can be far from the exact eigenvalues of A even if A is a floating-point matrix nearest to X −1 SX. In addition, the cost of symbolic or rational arithmetic is much expensive compared to that of numerical computations. We propose methods that produce a matrix with specified eigenvalues on the basis of an error-free transformation of matrix multiplication. For a matrix S ∈ F n×n given by a user, we approximate S by S ∈ F n×n , and obtain A := X −1 S X ∈ F n×n using the proposed method. We prove that no rounding error occurs in evaluating a product of three matrices X −1 S X using floating-point arithmetic. Hence, the exact eigenvalue of A can be known from the matrix S . This paper is organized as follows: In Section 2, basic theorems for floating-point numbers are introduced. In Section 3, we propose an error-free transformation for a product of three floating-point matrices. In Section 4, several methods for generating test matrices with exactly known eigenvalues are proposed based on the theory in Section 3, and we introduce numerical examples in checking the accuracy of numerical computations.

Basics of floating-point arithmetic
Notation and basic lemmas on floating-point arithmetic are first introduced. Notations fl(·), fl (·), and fl (·) indicate a computed result by floating-point arithmetic with rounding to nearest mode (roundTiesToEven), rounding downward mode (roundTowardNegative), and rounding upward mode (roundTowardPositive), respectively. Such rounding modes are defined in IEEE 754 [5]. We omit to write fl(·) for each operation for simplicity of description. For example, if we write fl((a + b) + c) for a, b, c ∈ F, then it means fl(fl(a + b) + c). Assume that neither overflow nor underflow occurs in fl(·) and fl (·). Let u and u S be a roundoff unit and the smallest positive number in F, respectively. A constant realmax is the largest floating-point number in F. For example, u = 2 −53 , u S = 2 −1074 , and realmax = 2 1024 (1 − u) for binary64 in IEEE 754. A function ufp(a) for a ∈ R, the unit in the first place of the binary representation of a, is defined as From the definition of the ufp function, we have We introduce several basic facts of floating-point numbers and floating-point arithmetic for a, b ∈ F and x ∈ R: ufp(a + b) ≤ ufp(fl(a + b)).
Applying Theorem 3.5 in [8] and an induction of binary tree described in [6], we can immediately obtain the following lemma.
|p i | ≤ σ , then no rounding error occurs in the sum by any orders.
We assume that divide and conquer methods by, for example, Strassen [9] and Winograd [11] are not applied for dot product, matrix-vector, and matrix-matrix multiplication.

General theory for error-free computations
For Y ∈ F m×n , S ∈ F n×p , X ∈ F p×q , we obtain a matrix S ∈ F n×p (S ≈ S) in order to compute Y S X without a rounding error. This technique is applied to generation of test matrices with exactly known eigenvalues.
We set several constants as follows: We define φ ij for x ij = 0 and ψ ij for y ij = 0 such that where all φ ij and ψ ij are powers of two. The constants φ ij and ψ ij imply the unit in the last non-zero bit of x ij and y ij , respectively. Exceptionally, we set φ ij = 0 for x ij = 0 and ψ ij = 0 for y ij = 0. Next, we set constants β, γ , θ, ω as follows: These imply for all (i, j ) pairs. Note that β, γ , θ, ω ≥ 1 and max both X and Y are non-singular matrices.
If s ij = 0 for all (i, j ) pairs, then it is trivial that the matrix A is the zero matrix. Hence, we assume max i,j |s ij | > 0. Let n Y : the maximum of the number of non-zero elements in row of Y , -n S : the maximum of the number of non-zero elements in row of S, -n X : the maximum of the number of non-zero elements in column of X.
We set n = min(n S , n X ). We define a constant σ as where Then, is satisfied because all of β, γ , θ, and ω are powers of two. It is possible to set α as fl (n · n · max i,j |s ii |) using floating-point arithmetic. A matrix S is computed by We introduce several lemmas. (14), β and γ in (11), and ψ and ω in (12), and are satisfied for all (i, j ) pairs.

Lemma 5 Assume that s ij is obtained by
Proof Using the definition of σ in (14) gives σ ≥ |s ij | for all (i, j ) pairs. From this and (7), we derive (6) and Lemma 4. Therefore, from Lemma 3 We show that rounding error never occurs in the evaluation of Y S X.
Hence, we assume there exists x ij = 0 for any j for the rest of the proof. Since n Y , n , β, θ, ω ≥ 1 and 4n Y · n · uβγ θω ≤ 1 from the assumption, we have Using Lemma 1 and (10) gives and they yield Next, we take an upper bound of p k=1 s ik x kj using Lemma 5, Lemma 3, and (24): Therefore, from (26) and (28), Lemma 2 derives fl(S X) = S X.
s k x kj .
If we change the definition of n , it is possible to verify fl((Y S )X) = Y S X similarly. Remark that the proposed method has reproducibility, 2 because there is no rounding error in any order of computations in fl(Y (S X)).

Test matrices based on specified matrices
In this section, we introduce how to generate a matrix with exactly known eigenvalues. We prepare S, X, Y ∈ F n×n (Y = X −1 ). Assume that exact eigenvalues can be known quickly from the matrix S. For example, S is a diagonal matrix, triangular matrices, and tridiagonal Toeplitz matrix. We use MATLAB R2020a on Windows 10 with Intel(R) Core(TM) i7-8665U in numerical examples in this section.
We introduce a basic algorithm for the generation of test matrices using MATLABlike code on the basis of the discussion in Section 3. Greek alphabets cannot be used for variables on MATLAB. However, we exceptionally accept to use Greek alphabets in the description of algorithms for consistency with the theories in this paper.
In the algorithm, ufp(·) at line 11 is not a built-in function in MATLAB. The calculation of ufp(a) for a ∈ F involves only three floating-point operations (see Algorithm 3.6 in [7]). Therefore, the cost for the ufp function is negligible. The function get bg at line 6 in Algorithm 1 produces β and γ according to (11). The function get tw similarly gives θ and ω based on (12).
To set multiple eigenvalues can be possible, because if s ij = s k , then z ij = z k . We can also generate a rank deficient matrix with exactly known eigenvalues. However, even if s ij = s k , z ij = z k may be satisfied. It happens for relatively small |s ij | compared to max i,j |s ij |. If uσ ≥ |s ij |, then z ij = 0. In the worst case, the matrix Z becomes the zero matrix, and an obtained matrix A is useless. In addition, a test matrix with the specified singular values can be generated using Y = X T .
Note that we can generate a real matrix with complex eigenvalues. Let i be the imaginary unit. The basic of linear algebra derives that a ± bi are eigenvalues of the following skew-symmetric matrix: If one wants to obtain a matrix with exactly known eigenvectors, then we obtain the matrix as follows.
in some cases due to a problem of rounding error. Then, S is not a skew-symmetric matrix. To solve the problem, we compute only Z is not a skew-symmetric matrix, while S is a skew-symmetric matrix.
The condition X, Y ∈ F n×n (Y = X −1 ) is very special. Therefore, we use specified matrices for X in (1) from the following sections.

Using weighing matrix
Assume that D ∈ F n×n is a diagonal matrix (n = 1). First, we use a weighing matrix. Let I be the identity matrix. If W ∈ R n×n with weight c ∈ N is a weighing matrix, W T W = cI and w ij ∈ {−1, 0, 1} for all (i, j ) pairs. The weighing matrix is useful for the candidate of the matrix X, because β = γ = θ = ω = 1 in (11) Then, we can compute W T n,c (D W n,c ) =: A without rounding errors. Hence, cd ii is the exact eigenvalue of the matrix A. Note that cd ii may not be representable in the floating-point numbers. An error-free transformation [3] produces floating-point vectors p and q such that It is well known that we quickly obtain p and q using fused multiply-add (FMA instruction defined in IEEE 754). The function fma (a, b, c) for a, b, c ∈ F provides a nearest floating-point number of ab + c. Then, Hadamard matrix, the special case of weighing matrix (c = n), is very useful for the candidate of X.
The following codes generate a test matrix on MATLAB. Let d ∈ F n be a vector with expected eigenvalues. The exact eigenvalue of A ∈ F n×n is p i + q i ≈ d i . Note that the function hadamard on MATLAB R2020a in the following algorithm works when n, n/12 or n/20 is a power of two. 3 Note that Algorithm 2 produces a real symmetric matrix. For n ∈ F and y ∈ F n , the function TwoProduct at line 12 transforms ny to p + q, p, q ∈ F n . If n is a power of two, q at line 12 is the zero vector.
At line 11, a matrix multiplication is computed. In the special cases, we can reduce the cost of the matrix multiplication. Let H n be a Hadamard matrix with dimension n. We can generate one of 2n-by-2n Hadamard matrices by the following manner [10] For a diagonal matrix D ∈ F 2n×2n , there exist matrices M, N ∈ R n×n such that where D 1 and D 2 are n-by-n diagonal matrices. Then, the matrices M and N can quickly be obtained, because D 1 and D 2 are diagonal matrices. Finally, we compute H T 2n DH 2n by Therefore, we only compute H T n M = H T n D 1 H T n and H T n N = H T n D 2 H T n as matrix multiplications to obtain H T 2n DH 2n . The cost of matrix multiplications can be reduced using this technique recursively.
Let r · n be a matrix size for the output, where r, n ∈ N and r is a power of two.
where D i ∈ F n×n , 1 ≤ i ≤ r are diagonal matrices, and O is the zero matrix. Then, we only compute the following r matrix multiplications: and C ij is obtained by The total cost becomes 2kn 3 + O(n 2 + k) floating-point operations for r-by-r block division. Note that the original cost without block computations is 2(kn) 3 +O((kn) 2 ), that is much expensive compared to block computations for large k.
We introduce an algorithm producing a matrix with eigenvalues based on the block computations.
We show numerical examples for the performance using this technique. Figures 1  and 2 show computing times with r-by-r block computations. r = 1 means a computing time for Algorithm 2. Computing times for solving linear systems (linsolve) and eigenvalue problems (eig) using MATLAB are shown for comparison. Note that eig in this example produces only approximate eigenvalues such that d = eig(A). The figures indicate that computing time for the generation of test matrices is much shorter than that for solving linear systems and that for eigenvalue problems.

Remark 1
We mentioned that MATLAB's hadamard function works in the case that n, n/12, or n/20 is the power of two. Even if n does not satisfy that restriction, the test matrices can be generated. Assume that H k exists and k < n. Let P ∈ R (n−k)×(n−k) be a permutation matrix. We generate instead of the Hadamard matrix. Then, D is obtained by The exact eigenvalue is kd ii for 1 ≤ i ≤ k and d ii for k + 1 ≤ i ≤ n.
We checked the accuracy of approximate eigenvalues obtained by eig. The vector d ∈ F n , n = 4096 in Algorithm 2 is set as follows. Ex.1: We compare accuracy of eigenvalues by d = eig(A) and [V , D] = eig(A) in MATLAB. Letλ 1 ≤λ 2 ≤ · · · ≤λ n be computed eigenvalues. Figures 3 and 4 show relative errors

Eigenvector matrix X using random integer matrices
We use the LU decomposition and generate integer matrices 4 for X. Let L and U be a unit lower triangular matrix and a unit upper triangular matrix, respectively. We generate X = LU , where we set L and U as 0-1 matrices in order to obtain small constants β, γ , θ, and ω. Then, Y = X −1 is also an integer matrix. Because L and U are 0-1 matrices, there is no rounding error in fl(LU ) up to n = u −1 . On the other hand, the problem is that the exact Y = X −1 can be obtained or not. LetỸ be a numerically obtained result for X −1 . I is the identity matrix, thenỸ = X −1 . Using MATLAB notation, L and U are given as L = spones(tril(sprandn(n, n, k))); L(1 : n + 1 : n * n) = 1; (32) U = spones(triu(sprandn(n, n, k))); U(1 : n + 1 : n * n) = 1; (33) where 0 ≤ k ≤ 1 is density of a matrix. The sparsity of L and U is controlled by k ∈ N.
The following algorithm generates A using random integer matrices. We introduce numerical examples for Algorithm 5. We set n = 1, 000, 5, 000, and d ∈ F n is generated by randn(n, 1) using MATLAB R2020a. Table 1 shows -The product of four constants βγ θω -The 2-norm condition number of X, i.e., X −1 2 X 2 = Y 2 X 2 -The relative change between d i and r i , i.e., max i |d i − r i | |d i | for each k. The numbers in Tables 1 and 2 are averages of 10 examples. If k is small, we can obtain small βγ θω. It makes the relative change for v small. If k becomes large, then -The relative change for a small absolute value, -The condition number of X, and -Density of A also become large.

Conclusion
We proposed methods generating matrices with specified eigenvalues based on the error-free transformation of floating-point numbers. They are useful for checking the accuracy and stability of numerical algorithms. Finally, we provided several possibilities for an eigenvector matrix and showed how to generate a real matrix with complex eigenvalues. The discussion can be extended to the generation of test matrices with complex eigenvalues and eigenvectors.