ABS-Based Direct Method for Solving Complex Systems of Linear Equations

Abstract

Efficient solution of linear systems of equations is one of the central topics of numerical computation. Linear systems with complex coefficients arise from various physics and quantum chemistry problems. In this paper, we propose a novel ABS-based algorithm, which is able to solve complex systems of linear equations. Theoretical analysis is given to highlight the basic features of our new algorithm. Four variants of our algorithm were also implemented and intensively tested on randomly generated full and sparse matrices and real-life problems. The results of numerical experiments reveal that our ABS-based algorithm is able to compute the solution with high accuracy. The performance of our algorithm was compared with a commercially available software, Matlab’s mldivide (\) algorithm. Our algorithm outperformed the Matlab algorithm in most cases in terms of computational accuracy. These results expand the practical usefulness of our algorithm.

1. Introduction

The problem of solving systems of linear equations plays a central role in diverse scientific fields such as signal processing, economics, computer science, and physics [1,2,3,4,5]. Often the problems arising from areas of practical life result in a system of equations with real coefficients, but there are important applications that lead to the following complex linear systems:

$$Ax=b,\mathrm{where}A\in {\mathbb{C}}^{m,n},x\in {\mathbb{C}}^n\mathrm{and}b\in {\mathbb{C}}^m.$$

(1)

Partial differential equations modelling dissipative processes usually involve complex coefficient functions or complex boundary conditions [6]. Other applications leading to complex linear systems include the discretization of time-dependent Schrödinger equations with implicit differential equations [7,8,9], inverse scattering problems, underwater acoustics, eddy current calculations [10], diffuse optical tomography [11], numerical calculations in quantum chromodynamics and numerical conformal mapping [12]. There are several methods to solve complex systems of linear equations. Moreover, a complex linear system with n unknowns can be reformulated to a linear system of equations with $2n$ real coefficients [13].

There are two basic popular approaches in solving complex linear systems of equations. The first one is when the preconditioned classical conjugate gradient method is used to solve the system of equations [14,15,16]. In most of these cases, the algorithms often do not work with the original general, i.e., non-Hermitian A coefficient matrix, but with a modified Hermitian positive definite normal equations

$$A^{H}Ax=A^{H}b,\mathrm{where}{A}^H\mathrm{is}\mathrm{the}\mathrm{conjugate}\mathrm{transpose}\mathrm{of}\mathrm{the}A\mathrm{matrix}.$$

(2)

The second popular approach is to solve usually non-symmetric linear systems of equations using one of the generalized CG methods [17,18], such as GMRES ([19] pp. 164–185) or other approaches based on the Arnoldi or the Lanczos biconjugate algorithms. These approaches generate an orthogonal basis for the Krylov subspaces associated with A and an initial vector $v_1$. The usual way to obtain an approximation to the exact solution of (1) from this basis is to force a biconjugate gradient condition. In both cases, the resulting iterative methods tend to converge relatively slowly.

Mainly, two factors influence the practical usability of a method: the accuracy of the solutions and the computation cost of the method. There are two main classes of techniques for solving linear systems of equations: iterative and direct methods.

Iterative solvers are generally faster on large-scale problems, while direct ones give more accurate results. Not surprisingly, iterative methods have come to the fore in solving complex linear systems of equations and many new algorithms have been published in recent decades [20,21,22]. However, the robustness and the speed of the convergence of iterative algorithms significantly depend on the condition number of the coefficient matrices. To avoid convergence problems, robust preconditioners need to be used. In some cases, problem-specific preconditioners are highly effective, but they are often difficult to parallelize on modern high performance computing (HPC) platforms [23]. Moreover, Koric et al. [24] revealed that no iterative preconditioned solver combination could correctly solve the highly ill-conditioned systems of equations. Therefore, researchers often have to compromise between accuracy of the solutions and computational effort. Direct solvers typically do not have these limitations but require much more computational power to execute. Thus, today’s high performance solutions such as parallel computing frameworks, may provide new possibilities for direct algorithms and they can be more suitable choices to ill-conditioned problems.

In this paper, we present a new ABS-based direct method that can solve large-scale problems. The ABS class of methods was originally introduced by Abaffy, Broyden and Spedicato (1984) developed to solve linear systems of equations where coefficients of the equations are real numbers [25]. These algorithms can also be used for various other purposes such as solving non-linear systems of equations or optimization problems [26,27,28]. ABS-based algorithms can be easily and effectively parallelised [29,30] which underlines the practical usefulness of these algorithms.

The remainder of this paper is organized as follows. In Section 2, we present our new scaled complex ABS algorithm, which can solve complex linear systems of equations and we prove some of its basic features. We also show a special choice of parameters of our algorithm, which ensures that the search vectors ($p_i$) are $A^{H}A$ conjugate. We provide a detailed numerical analysis of four variants of our algorithm in Section 3. Section 4 contains concluding remarks and a brief description of our plans on the application of the outcomes of this work.

2. The Scaled Complex ABS (scABS) Algorithm

In this section, we present our new ABS-based algorithm (scABS) that can solve systems of complex linear problems. Instead of the original system (1) let us consider the following scaled complex system of linear equations

$$V^{H}Ax=V^{H}b,\mathrm{where}V\in {\mathbb{C}}^{m,n}\mathrm{is}\mathrm{an}\mathrm{arbitrary},\mathrm{non}\text{-}\sin \mathrm{gular}\mathrm{matrix}.$$

(3)

Note that the systems (1) and (3) are equivalent. Thus, if we want our algorithm to solve the original unscaled system of Equation (1), then by choosing the scaling matrix as the unit matrix we obtain the desired formulas and values. The sole role of the scaling matrix V is to provide significant flexibility to the algorithm [25] which allows us to reduce the computational inaccuracy [31] and also to reduce the computational complexity of the algorithm [27].

It is also of interest to note that our scaled complex ABS algorithm defines not one single algorithm, but a class of variants. Each particular variant is determined by the choice of the parameters $H_1$, $v_i$, $z_i$, $w_i$. These variants have different properties, e.g., by choosing $H_1=I$ unit matrix, $v_i=A{H}_i^H{z}_i$, and $w_i$ to be a multiple of $z_i$, the ABS algorithm is a reformulation of the QR algorithm, as we will discuss in Section 2.2.

The state of the complex system of linear equations is checked by the variable $s_i=H_i{A}^H{v}_i$. An important property of our algorithm is that $s_i$ is zero if, and only if, the current row of matrix A is a linear combination of the previous rows. It depends on the value of the right-hand side whether our system is then a linear combination of the previous equations or even incompatible with them. We use the value of $x_i$ to distinguish between the two cases, namely whether it solves the $ith$ equation or not. If it does, we simply skip the equation ($x_{i+1}=x_i$, $H_{i+1}=H_i$), otherwise, we stop the algorithm. The state of the system is stored in the $iflag$ variable. If $iflag$ is zero, then neither linear dependency nor incompatibility is detected in the algorithm. If the value of $iflag$ is positive, then a number of linearly dependent equations are found. If the value of $iflag$ is negative then the $-iflagth$ equation is incompatible with the previous ones.

$H_1$ was selected to be a unit matrix for the sake of simplicity. However, $H_1$ could theoretically be any arbitrary unimodular nonsingular matrix provided that $H_1\in {\mathbb{C}}^{n,n}$.

2.1. Basic Properties of the Scaled Complex ABS Algorithm

We consider some fundamental properties of the $H_i$ generated by the scaled complex ABS (scABS) algorithm (Algorithm 1).

Theorem 1.

Given the $A^H{v}_1$,..., $A^H{v}_m$ $\in {\mathbb{C}}^n$ vectors and the $H_1\in {\mathbb{C}}^{n,n}$ arbitrary nonsingular matrix, consider the sequence of matrices $H_2,\dots ,H_{m+1}$ generated by (6). The following relations are true for $i=2,\dots ,m+1$:

$$H_i{A}^H{v}_j=0,1\le j<i.$$

(4)

Proof. 

We proceed by induction. For $i=2$, the theorem is true since

$$H_2=H_1-\frac{H_1{A}^H{v}_1\ast w_1^H{H}_1}{\widetilde{w_1}}\overline{w_1^H{H}_1{A}^H{v}_1}.$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_2{A}^H{v}_1& =H_1{A}^H{v}_1-\frac{H_1{A}^H{v}_1\ast w_1^H{H}_1}{\widetilde{w_1}}\overline{w_1^H{H}_1{A}^H{v}_1}\ast A^H{v}_1\hfill \\ & =H_1{A}^H{v}_1-\frac{H_1{A}^H{v}_1\ast w_1^H{H}_1{A}^H{v}_1}{\widetilde{w_1}}\overline{w_1^H{H}_1{A}^H{v}_1}\hfill \\ & =H_1{A}^H{v}_1-\frac{H_1{A}^H{v}_1\ast \widetilde{w_1}}{\widetilde{w_1}}=0.\hfill \end{array}\end{array}$$

Assuming that the theorem is true up to $i<m$, we prove it for $i+1$.

$$H_{i+1}=H_i-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_{i+1}{A}^H{v}_j& =H_i{A}^H{v}_j-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\ast A^H{v}_j\hfill \\ & =H_i{A}^H{v}_j-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i{A}^H{v}_j}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & =H_i{A}^H{v}_j-\frac{H_i{A}^H{v}_i\ast \widetilde{w_i}}{\widetilde{w_i}}=0.\hfill \end{array}\end{array}$$

   □

Algorithm 1: Scaled complex ABS (scABS) algorithm.

Remark 1.

The null space of complex projection matrices $Null\left(H_i\right)=\left\{A^H{v}_1,A^H{v}_2,\dots ,A^H{v}_{i-1}\right\}$ follows immediately from Theorem 1 and it means that $H_{m+1}=0$

Remark 2.

$H_i{A}^H{v}_i$ computed by the scABS algorithm is zero only if $A^H{v}_i$ is a linear combination of $A^H{v}_1,A^H{v}_2,\dots ,A^H{v}_{i-1}$.

Remark 3.

The scaled complex ABS algorithm is well defined.

Theorem 2.

Consider the matrices $H_i$ generated by (6) with the starting $H_1$ is nonsingular. For $i=1,\dots ,m+1$ and $1\le j\le i$ the following relations are true:

$$H_i{H}_1^{-1}{H}_j=H_i,$$

(7)

$$H_j{H}_1^{-1}{H}_i=H_i.$$

(8)

Proof. 

We only prove (7) since the proof for (8) is similar. We proceed by induction. For $i=1$, $H_1{H}_1^{-1}{H}_1=H_1$ is trivially true. Assuming now that (7) is true up to the index i, we have to prove that $H_{i+1}{H}_1^{-1}{H}_j=H_{i+1}$.

For $j=i+1$, we have

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill H_{i+1}{H}_1^{-1}{H}_{i+1}& =(H_i-\frac{H_i{A}^H{v}_1\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i})H_1^{-1}\hfill \\ & \ast (H_i-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i})\hfill \\ & =H_i{H}_1^{-1}{H}_i-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i{H}_1^{-1}{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & -\frac{H_i{H}_1^{-1}{H}_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & +\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}{H}_1^{-1}\ast \frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & =H_i{H}_1^{-1}{H}_i-2\ast \frac{H_i{A}^H{v}_i\ast w_i^H{H}_i{H}_1^{-1}{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}+\hfill \\ & \frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}{H}_1^{-1}\ast \frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}{\left(\overline{w_i^H{H}_i{A}^H{v}_i}\right)}^2\hfill \\ & =H_i-2\ast \frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & +\frac{H_i{A}^H{v}_i\ast \widetilde{w_i}\ast w_i^H{H}_i}{{\widetilde{w_i}}^2}\overline{w_i^H{H}_i{A}^H{v}_i}\hfill \\ & =H_i-\frac{H_i{A}^H{v}_i\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}=H_{i+1}.\hfill \end{array}\end{array}$$

For $j<i+1$, we have

$$\begin{array}{cc}\hfill H_{i+1}{H}_1^{-1}{H}_j& =(H_i-\frac{H_i{A}^H{v}_1\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i})H_1^{-1}{H}_j\hfill \\ & =H_i-\frac{H_i{A}^H{v}_1\ast w_i^H{H}_i}{\widetilde{w_i}}\overline{w_i^H{H}_i{A}^H{v}_i}=H_{i+1}\hfill \end{array}$$

again by the induction.    □

Remark 4.

If $H_1$ is the unit matrix then Theorem 2 implies that $H_i$ matrices are idempotent, i.e., $H_i{H}_i=H_i$.

Remark 5.

The $H_i{H}_i=H_i$ feature of our algorithm enables further potential improvement in computational accuracy of the scaled complex ABS algorithm. We have recently published that the reprojection of real ABS algorithms, i.e., using $p_i=H_i^H(H_i^H{z}_i$) projection vectors instead of $p_i=H_i^H{z}_i$, enhances the precision of numerical calculation [30,31], which expands the practical usefulness of our algorithms.

Theorem 3.

Any vector $y\in {\mathbb{C}}^n$ that solves the first i equations of (3) can be formulated as

$$y=x_{i+1}+H_{i+1}^{H}s$$

(9)

where s is a vector in ${\mathbb{C}}^n$.

Proof. 

Decompose $y\in {\mathbb{C}}^n$ as $y=x_{i+1}+y^{+}$. As y solves the first i equations,

$$A^H{v}_{j}y=b^H{v}_j=A^H{v}_j{x}_{i+1}+A^H{v}_j{y}^{+},\mathrm{where}1\le j\le i.$$

(10)

Thus, $A^H{v}_j{y}^{+}=0$ for $1\le j\le i$ which means that $y^{+}\perp \{A^H{v}_1,..,A^H{v}_i\}=\perp Null$ $\left(H_{i+1}\right)$, so $y^{+}=H_{i+1}^{H}r$ where $r\in {\mathbb{C}}^n$ according to the properties of the scaled complex ABS algorithm. So $y^{+}=H_{i+1}^{H}r$ where $r\in {\mathbb{C}}^n$ according to the properties of ABS.

$$y^{+}=H_{i+1}^{H}r=H_{i+1}^H·H_{i+1}^{H}r=H_{i+1}^H{y}^{+}\mathrm{where}{y}^{+}\in {\mathbb{C}}^n.$$

(11)

   □

Remark 6.

Referring to the proofs of the original ABS algorithm [26], it can be shown that, for simplicity, assuming that the A matrix is full rank, $L_i=V_i^H{A}_i{P}_i$ is a nonsingular lower triangular matrix where $P_i=(p_1,\dots ,p_i)$, $V_i=(v_1,\dots ,v_i)$ and $A_i=(a_1^H,\dots ,a_i^H)$ the transpose of the rows of the matrix A. If $i=m$, then the following semifactorization of the inverse can be obtained:

$$A^{-1}=P{L}^{-1}{V}^H,\mathit{where}P=P_m,\mathit{and}V=V_m.$$

(12)

This semi-factorization of the inverse of the A matrix may provide an opportunity to develop an ABS-based preconditioner to accelerate various Krylov subspace methods. For several choices of the matrix V, the matrix L is diagonal, hence formula (12) gives an explicit factorization of the inverse of the coefficient matrix A.

Remark 7.

Examining the properties of the scaled complex ABS algorithm in terms of practical usability, we have already mentioned in Remark 4 that the idempotent property of the $H_i$ matrix can be used to increase the computational accuracy of our algorithm by reprojection stepwise. Broyden showed [32] that the computation error in the solution is reduced up to 2 orders in $\frac{\parallel A^T{v}_i{\parallel }_2}{\parallel H_i{A}^T{v}_i{\parallel }_2}$ for the real ABS algorithm. If the projection vectors ($p_i$) are re-projected with an additional computational cost, i.e., $p_i=H_i^T{z}_i=H_i^T{H}_i^T{z}_i$, a more accurate result can be obtained due to the cancellation errors in computational accuracy [33,34]. Interestingly, our preliminary results showed that the scABS algorithm has similar properties. Increasing accuracy in this way will result in a significant increase in computational costs. One solution to speed up the algorithm may be to parallelize the processes. In a recent paper [30], we found that for a real ABS algorithm, parallelization yielded a significant gain in computational speed.

2.2. Special Choice of the Scaling Vector

We consider the following choice of the scaling vector:

$$v_i=A{p}_i.$$

(13)

This selection of the scaling vector has many beneficial features such as

• the scaling vectors defined by (13) are mutually orthogonal, i.e.,
  $v_i^H{v}_j=0$ for $i>j$ and $v_i^H{v}_i\ne 0$.
• the search vectors $p_i$ are $A^{H}A$ conjugate if A is square and nonsingular matrix.
• the algorithm generates an implicit factorization of A into the product of an orthogonal and an upper triangular matrix if $w_i$ is a multiple of $z_i$ and $H_1=I$ (unit matrix). The algorithm with these selections can be considered as an alternative method of the classic QR factorization.
• by the choice of $H_1=I$, $z_i=w_i=e_i$, the required number of multiplication is $\frac{11}{6}{n}^3+\mathcal{O}\left(n^2\right)$.

The above statements about of the oscABS algorithm (Algorithm 2) can be easily verified based on the proofs of Theorem 8.11, Theorem 8.16 and Theorem 8.18 in [25].

Algorithm 2: Orthogonally scaled complex ABS (oscABS) algorithm.

3. Numerical Experiments

We were also interested in the numerical features of our scaled complex ABS algorithm. To this end, four variants of the orthogonally scaled complex ABS algorithm were implemented in MatlabR2013b (Mathworks, Inc., USA). The experiments were performed on a personal computer with Intel Core i7-2600 3.4GHz CPU with integrated graphics, 4 GB RAM running Microsoft Windows 10 Professional and MATLAB version R2013b. No software other than the operating system tasks, MATLAB and ESET NOD32 antivirus were running during the experiments.

The four implemented variants of the oscABS algorithm and the Matlab function used are as follows:

• S3rr: $z_i$ and $w_i$ are selected such that $z_i=r_i$, $w_i=r_i$.
• S3ATA: $z_i$ and $w_i$ are selected such that $z_i=A^H{r}_i$, $w_i=A^H{r}_i$.
• S3ee: $z_i$ and $w_i$ are selected such that $z_i=e_i$, $w_i=e_i$.
• S3ep: $z_i$ and $w_i$ are selected such that $z_i=e_i$, $w_i=H_i^H{p}_i$.
• Matlab: the built-in mldivide (\) function of Matlab.

We evaluated all methods in aspects of accuracy of the computed solution $\parallel A{x}_n{-b\parallel }_2$, and execution time in seconds.

3.1. Randomly Generated Problems

In our first numerical experiments, we tested the variants of our algorithm on randomly generated dense and sparse matrices. These matrices were generated using the rand and sprand MATLAB functions. In general, we performed 10 separate tests on independent matrices to obtain each data point in this subsection. The mean values from those results are depicted in the following figures. The standard variation of the values usually fell within well below 5% of the mean value (data not shown). The solutions were randomly generated using the rand function and the right sides of the systems were calculated as the products of the matrix and the solution. In these experiments, we tested how the performance of our algorithm was affected by gradually increasing the dimension of the coefficient matrix from $10\times 10$ to $1500\times 1500$. As shown in Figure 1, each of the variants of the ABS-based algorithm was able to solve the systems of equations within ${10}^{-10}$ accuracy, and overall, the S3ee variant was the most accurate. We did not see a large difference up to 1500 dimensions between the four variants of the orthogonally scaled complex ABS algorithm.

Our next aim was to compare the performance of the best variant of oscABS algorithm (S3ee) with a commercially available software able to solve complex linear systems of equations. We chose the mldivide (\) algorithm of the Matlab software package from MathWorks, Inc. for that purpose. Note that Matlab’s mldivide algorithm is not one specific algorithm, but consists of several algorithms, from which it chooses depending on the property of the matrix A. The mldivide function performs several checks on the coefficient matrix to determine whether it has some special property (e.g., whether it is sparse or symmetric) and, knowing this, selects the appropriate matching solver. For example, if A is full but not a square matrix, then it uses the QR algorithm, if not, then depending on the properties of A, possibly the Hessenberg, Cholesky, LU, or even the LDL solver. This means that we compared our algorithm with different solvers, selected to be most appropriate for the given problem by the Matlab program [35].

The experiments were performed with the randomly generated matrices described above.

Our results are summarized in Figure 2. The S3ee algorithm outperformed Matlab algorithm in both full rank, indefinit and rank-deficient problems as well. It is clear from Figure 2 that the difference increases significantly as the dimension increases.

Our next set of experiments focused on the analysis of the computational accuracy of the algorithms under different matrix densities. These experiments were performed with uniformly distributed random numbers of the coefficient matrix at different densities. As shown in Figure 3, the S3ATA implementation calculated the solution most accurately and the S3rr the least accurately in most cases. The other two implementations of the orthogonally scaled complex ABS algorithms worked similar to S3ATA (i.e., within the same order of magnitude). Even at 1500 dimensions, these algorithms were able to calculate the result with an accuracy of ${10}^{-11}$. It is also worth noting that reducing the density from $100\%$ to $1\%$ resulted in a significant improvement in computational accuracy. The accuracy increased from ${10}^{-11}$ to ${10}^{-13}$ for the S3ATA, S3ee, and S3ep variants.

Next, we compared the computational accuracy of our ABS-based algorithm with the Matlab algorithm on problems of different densities. As shown in Figure 4, the S3ee variant of orthogonally scaled complex ABS algorithms outperformed the Matlab algorithm again. It is worth noting that in the case of the Matlab algorithm, the change in matrix density resulted in a moderate improvement in the accuracy of the computations, as it increased from ${10}^{-11}$ to ${10}^{-12}$.

In order to gain a deeper understanding of the numerical characteristics of the four variants of our algorithm, we compared the computational accuracy and execution time required to solve the 1500-dimensional systems of equations (Figure 5).

Figure 5 shows that the numerical properties of the S3ee variant are the best, both in terms of computational accuracy and execution time. This may be explained by the fact that no computation is needed for parameter definitions ($z_i,w_i$) since the appropriate unit vectors are used and the lack of computation may also explain the accuracy since round-off errors do not accumulate. This round-off error may also explain the relatively poor numerical performance of S3ATA and S3rr.

To obtain a more comprehensive view of the numerical properties of our algorithms, we also investigated their computational accuracy for large, i.e., more than 5000 dimension problems. The dimensions of these problems were determined according to Golub and Van Loan’s suggestion [36] that for each problem, the q value, which is the characteristic of the problem, remains substantially below 1. The $q=u·k_{\infty }\left(A\right)$, where u (The value of the unit round-off in Matlab is $2^{-52}$ ) is unit round-off and $k_{\infty }\left(A\right)$ condition number of the coefficient matrix A in infinity-norm. Our results for the high-dimensional problems are summarized in Figures S1 and S2 in the Supplementary Material. Our experiments in high dimensions clearly showed that the different variants of the ABS-based algorithm solved the random problems significantly more accurately than Matlab solver. This is especially true in the rank-deficient cases, where for 5000 dimensions the Matlab function was not able to solve any problem.

3.2. Chosen Test Problems from MATLAB Gallery

Our next experiments compared the computation accuracy of our algorithm on known complex matrices in the Matlab Gallery.

Table 1. Summary of test problems chosen from MATLAB Gallery.

Name	Description
Symmetric Toeplitz	Diagonal-constant matrix
	using Matlab toeplitz(r) function where r is the first row.
	The vector r was randomly generated in this experiment.
Non-symmetric Toeplitz	Diagonal-constant matrix
	using Matlab toeplitz(c,r) function where c is first column,
	r is the first row of the matrix.
	The vector c, and the vector r were generated randomly.
Hankel	Symmetric and constant across the anti-diagonals matrix
	using Matlab hankel(c) function where c
	defines the first column of the matrix.
	The c column was generated randomly in this experiment.
Smoke	Complex matrix with a ‘smoke ring’ pseudospectrum
	using Matlab gallery(‘smoke’,n) function
	where n is the dimension of the matrix.

summarizes the problems we chose.

For each matrix, we tested the dimensions between 10–1500. Experiments presented in Figure 6 revealed that all four variants of the complex ABS algorithm were able to solve the problems with an accuracy of approximately ${10}^{-10}$ and in most cases the computation error remained significantly below ${10}^{-11}$. Examining the different variants of the orthogonally scaled complex ABS algorithm, we found that the S3ee and S3ep algorithms computed the problems with similar accuracy, with the S3ATA variant being the least accurate except for the Smoke problem.

Next, we compared the S3ee variant with the Matlab algorithm (see Figure 7). The ABS-based algorithm slightly outperformed the Matlab algorithm on all but the Smoke matrices.

In order to gain a deeper understanding of the numerical characteristics of the four variants of our algorithm, we compared the computational accuracy and execution time required to solve the 1500-dimensional systems of equations (see Figure 8).

For Matlab Gallery problems, the running properties of the algorithms are slightly different from those seen for randomly generated problems. For the first three matrix (Symmetric Toeplitz, Non-symmetric Toeplitz, Hankel) problems, we see similar results as for the randomly generated problems, the ABS-based variants giving more accurate results than the Matlab algorithm. However, for the fourth, Smoke matrix, for the first time Matlab’s algorithm gave the most accurate result tested. This fact can be partly explained by the special structure of the Smoke matrix since the diagonal of the n-dimensional Smoke matrix consists of the set of all nth roots of unity, and 1’s on its superdiagonal and 1 in the (n,1) position. It is interesting to note that the Matlab algorithm and the S3ATA variant behaved very similarly in these problems, both computing the solution with relatively large errors for the first three problems, but giving the most accurate solutions for the Smoke matrices. In the case of running speeds, it is clear that the lower computational demand of the S3ee variant in the 1500 dimension already resulted in significantly shorter running times.

We have also examined the behaviour of our ABS-based variants on large, i.e., more than 5000 dimensions of selected Matlab Gallery problems. Our results of the ABS-based variants are summarized in Figure S3, while the comparison of S3ee and Matlab results are summarized in Figure S4 in Supplementary Materials. For high dimensions, the numerical properties of the algorithms were very similar to those for low dimensions. Overall, of the ABS-based variants, the S3ee variant computed the most accurately and all ABS-based variants gave more accurate results than the Matlab algorithm except for the Smoke problem.

3.3. Real-Life Problems

We also wanted to examine how the four variants of our ABS-based algorithm work on real-life problems. Our next experiments focused on the performance of the algorithms on examples from the SuiteSparse matrix collection [37] and from the SWUFE-Math Test Matrices Library [38].

Table 2. Summary of used real-life problems. Matrices marked by italic font come from the SWUFE-Math Test Matrices Library, while non-italicized matrices come from the SuiteSparse matrix collection.

Matrix Name	n	m	Condition Number	Application Area
qc324	324	324	7.38 × ${10}^4$	Electromagnetics
young1c	841	841	9.91 × ${10}^2$	Acoustics
young2c	841	841	9.91 × ${10}^2$	Duplicate Acoustics
young4c	841	841	5.55 × ${10}^2$	Acoustics
dwg961b	961	961	3.35 × ${10}^7$	Electromagnetics
cube1800_test3	1800	1800	4.98 × ${10}^1$	Electromagnetics
parallelepipede_test2	2016	2016	8.43 × ${10}^1$	Electromagnetics
sphere2430_test1	2431	2431	4.32 × ${10}^1$	Electromagnetics
qc2534	2534	2534	5.19 × ${10}^5$	Electromagnetics
conf5_0-4x4-14	3072	3072	9.17 × ${10}^3$	Quantum Chemistry
mplate	5962	5962	4.83 × ${10}^{16}$	Acoustics

summarizes the matrices we used and their key properties.

The calculation accuracies are outlined in Panel A of Figure 9. It can be stated that each of the variants of the ABS-based algorithm was able to solve the problems with acceptable accuracy. For these problems, the S3rr algorithm performed best. It should be noted, however, that in most cases, the Matlab function calculated the solution most accurately. This may partly be explained by the fact that Matlab uses different solvers for different (i.e., sparse, dense) matrices. In addition to the accuracy of the solutions and execution times, we show the relative 2-norm of the residual vectors in the Panel C to ensure that the numerical properties of our ABS-based algorithm can be compared with other published algorithms [38,39,40] developed to solve linear systems of equations. These comparisons revealed that our method is significantly more accurate than expected for iterative solutions. Furthermore, some iterative methods achieved this accuracy of ${10}^{-6}$ only after a relatively large number of iterations. In the case of the young1c problem, several algorithms needed nearly 400 steps to achieve an accuracy of ${10}^{-6}$ [38,39], while our algorithm achieved an accuracy of ${10}^{-14}$ in about twice as many steps.

4. Conclusions

In this paper, we presented a new ABS-based algorithm for solving complex linear systems of equations and we proved some of its basic features. We also showed a special choice of parameters of our algorithm, which ensures that the search vectors ($p_i$) are $A^{H}A$ conjugate. A detailed numerical analysis of four variants of the ABS-based orthogonally scaled algorithm has also been provided. These variants were tested on randomly generated full and rank deficient, different-density systems of linear equations. Furthermore, the computational accuracy of the algorithm was tested on real-life problems. These numerical experiments showed that the ABS-based algorithm solved the problem with acceptable accuracy in all cases and provided more accurate results than the MATLAB built-in function in most test cases. These numerical results demonstrate the practical usefulness of the algorithm in addition to its theoretical significance.

Additionally, a valuable numerical property of our algorithm is that if we want to compute a complex system of linear equations with several right-hand sides [40], it is not necessary to recompute the matrix $H_i$. Instead, it is sufficient to store the vectors $p_i$ and recompute the updates $x_i$, which can significantly reduce the computational cost. and thus our algorithm can be used effectively to solve such problems.