Matrix LU Research Articles

Exact Matrix Factorization Updates for Nonlinear Programming

LU and Cholesky matrix factorization algorithms are core subroutines used to solve systems of linear equations (SLEs) encountered when solving an optimization problem. Standard floating-point algorithms are highly efficient but remain susceptible to the accumulation of round-off errors, which can lead solvers to return feasibility and optimality claims that are actually invalid. This paper introduces a novel direct solution approach for solving sequences of closely related SLEs encountered in nonlinear programming efficiently and without round-off errors. Specifically, it introduces rank-one update algorithms for the round-off error–free factorization framework, a tool set built on integer-preserving arithmetic that has led to the development and implementation of extremely reliable subroutines for solving SLEs occurring in linear programming. The formal guarantees of the presented algorithms are established through the derivation of theoretical insights. Their advantages are supported with computational experiments, which demonstrate upward of 75× improvements over exact factorization runtimes on fully dense matrices with more than one million entries. A significant advantage of the featured integer-preserving framework is that the length of any matrix coefficient produced by its algorithms is bounded polynomially in the size of the inputs without having to resort to greatest common divisor operations, which are required by and thereby hinder an efficient implementation of exact rational arithmetic approaches. History: Accepted by Antonio Frangioni, Area Editor for Design & Analysis of Algorithms–Continuous. Supplemental Material: The software that supports the findings of this study is available within the paper and its Supplemental Information ( https://pubsonline.informs.org/doi/suppl/10.1287/ijoc.2021.0331 ) as well as from the IJOC GitHub software repository ( https://github.com/INFORMSJoC/2021.0331 ). The complete IJOC Software and Data Repository is available at https://informsjoc.github.io/ .

Dynamical homotopy transient-based technique to improve the convergence of ill-posed power flow problem

This paper proposes a hybrid technique to solve the ill-posed Power Flow Problem (PFP), considering a homotopy approach. The primary proposal is to solve large-scale problems where the traditional Newton–Raphson (NR) fails to converge, as in the case of ill-posed systems. The method explores a dynamical homotopy transient-based technique to improve the convergence of the ill-conditioned problem instead of using the classical static method. Depending on the integration selected scheme and the integration step, the result furnished by the dynamical homotopy method has low accuracy. Then, the NR method is employed to refine the low-accuracy result and accurately determine the ill-posed PFP solution. The proposed approach can be implemented efficiently using only one Jacobian matrix computation and LU factorization per point of the homotopy path. In the static homotopy problem, a PFP using previous results must be solved per path point. In this case, some LU factorizations are necessary for each path point. The technique’s performance was evaluated through experiments, including a 70,000-bus large-scale system. The approximate dynamical homotopy result used as an initial estimate provided appropriate convergence quality for the NR method to determine a high-precision solution to the PFP.

Performance evaluation of LU matrix decomposition using the SYCL standard

The object of this study is the performance of the SYCL standard tools when solving the LU matrix decomposition problem. SYCL is a fairly new technology for parallel computing in heterogeneous systems, so the topic of evaluating the performance of the standard on specific tasks in the field of parallel computing is relevant. In the study, the algorithm of parallelized LU decomposition of a square matrix was implemented by means of the SYCL standard and standard C++, and an experiment was conducted to test the implementation in a heterogeneous system with several types of processors. During testing, the program received square matrices of various dimensions as input, and the output was the execution time of the LU schedule on the selected processor. The obtained results, presented in the form of tabular and graphic data, show the advantage of the implementation of the SYCL standard over ordinary C++ by more than 2 times when using a graphics processor. It was experimentally shown that the implementation on SYCL is almost not inferior in speed to the implementation on ordinary C++ when executed on a central processor. Such results are caused both by the high possibility of parallelizing the LU schedule algorithm itself, and by the great work of the developers of the standard on its optimization. The obtained results indicate the possibility of speeding up the solution of the LU decomposition of the matrix and similar algorithms by means of SYCL when using heterogeneous systems with processors optimized for data parallelism. The results of the study can be used in justifying the choice of technology for solving LU matrix decomposition problems or problems with a similar parallelization scheme.

Refined isogeometric analysis of quadratic eigenvalue problems

Certain applications that analyze damping effects require the solution of quadratic eigenvalue problems (QEPs). We use refined isogeometric analysis (rIGA) to solve quadratic eigenproblems. rIGA discretization, while conserving desirable properties of maximum-continuity isogeometric analysis (IGA), reduces the interconnection between degrees of freedom by adding low-continuity basis functions. This connectivity reduction in rIGA’s algebraic system results in faster matrix LU factorizations when using multifrontal direct solvers. We compare computational costs of rIGA versus those of IGA when employing Krylov eigensolvers to solve quadratic eigenproblems arising in 2D vector-valued multifield problems. For large problem sizes, the eigencomputation cost is governed by the cost of LU factorization, followed by costs of several matrix–vector and vector–vector multiplications, which correspond to Krylov projections. We minimize the computational cost by introducing C0 and C1 separators at specific element interfaces for our rIGA generalizations of the curl-conforming Nédélec and divergence-conforming Raviart–Thomas finite elements. Let p be the polynomial degree of basis functions; the LU factorization is up to O((p−1)2) times faster when using rIGA compared to IGA in the asymptotic regime. Thus, rIGA theoretically improves the total eigencomputation cost by O((p−1)2) for sufficiently large problem sizes. Yet, in practical cases of moderate-size eigenproblems, the improvement rate deteriorates as the number of computed eigenvalues increases because of multiple matrix–vector and vector–vector operations. Our numerical tests show that rIGA accelerates the solution of quadratic eigensystems by O(p−1) for moderately sized problems when we seek to compute a reasonable number of eigenvalues.

Performance impact of precision reduction in sparse linear systems solvers.

It is well established that reduced precision arithmetic can be exploited to accelerate the solution of dense linear systems. Typical examples are mixed precision algorithms that reduce the execution time and the energy consumption of parallel solvers for dense linear systems by factorizing a matrix at a precision lower than the working precision. Much less is known about the efficiency of reduced precision in parallel solvers for sparse linear systems, and existing work focuses on single core experiments. We evaluate the benefits of using single precision arithmetic in solving a double precision sparse linear system using multiple cores. We consider both direct methods and iterative methods and we focus on using single precision for the key components of LU factorization and matrix–vector products. Our results show that the anticipated speedup of 2 over a double precision LU factorization is obtained only for the very largest of our test problems. We point out two key factors underlying the poor speedup. First, we find that single precision sparse LU factorization is prone to a severe loss of performance due to the intrusion of subnormal numbers. We identify a mechanism that allows cascading fill-ins to generate subnormal numbers and show that automatically flushing subnormals to zero avoids the performance penalties. The second factor is the lack of parallelism in the analysis and reordering phases of the solvers and the absence of floating-point arithmetic in these phases. For iterative solvers, we find that for the majority of the matrices computing or applying incomplete factorization preconditioners in single precision provides at best modest performance benefits compared with the use of double precision. We also find that using single precision for the matrix–vector product kernels provides an average speedup of 1.5 over double precision kernels. In both cases some form of refinement is needed to raise the single precision results to double precision accuracy, which will reduce performance gains.

Statistical enumeration of groups by double cosets

Let H and K be subgroups of a finite group G. Pick g∈G uniformly at random. We study the distribution induced on double cosets. Three examples are treated in detail: 1) H=K= the Borel subgroup in GLn(Fq). This leads to new theorems for Mallows measure on permutations and new insights into the LU matrix factorization. 2) The double cosets of the hyperoctahedral group inside S2n, which leads to new applications of the Ewens's sampling formula of mathematical genetics. 3) Finally, if H and K are parabolic subgroups of Sn, the double cosets are ‘contingency tables’, studied by statisticians for the past 100 years.

Weighted Average Based Differential Quadrature Method for One-Dimensional Homogeneous First Order Nonlinear Parabolic Partial Differential Equation

In this paper, the weighted average-based differential quadrature method is presented for solving one-dimensional homogeneous first-order non-linear parabolic partial differential equation. First, the given solution domain is discretized by using uniform discretization grid point. Next, by using Taylor series expansion we obtain central finite difference discretization of the partial differential equation involving with temporal variable associated with weighted average of partial derivative concerning spatial variable. From this, we obtain the system of nonlinear ordinary differential equations and it is linearized by using the quasilinearization method. Then by using the polynomial-based differential quadrature method for approximating derivative involving with spatial variable at specified grid point, we obtain the system of linear equation. Then they obtained linear system equation is solved by using the LU matrix decomposition method. To validate the applicability of the proposed method, two model examples are considered and solved at each specific grid point on its solution domain. The stability and convergent analysis of the present method is worked by supported the theoretical and mathematical statements and the accuracy of the solution is obtained. The accuracy of the present method has been shown in the sense of root mean square error norm and maximum absolute error norm and the local behavior of the solution is captured exactly. Numerical versus exact solutions and behavior of maximum absolute error between them have been presented in terms of graphs and the corresponding root mean square error norm and maximum absolute error norm presented in tables. The present method approximates the exact solution very well and it is quite efficient and practically well suited for solving the non-linear parabolic equation. The numerical result presented in tables and graphs indicates that the approximate solution is in good agreement with the exact solution.

Weighted Average Based Differential Quadrature Method for One-Dimensional Homogeneous First Order Nonlinear Parabolic Partial Differential Equation

In this paper, the weighted average-based differential quadrature method is presented for solving one-dimensional homogeneous first-order non-linear parabolic partial differential equation. First, the given solution domain is discretized by using uniform discretization grid point. Next, by using Taylor series expansion we obtain central finite difference discretization of the partial differential equation involving with temporal variable associated with weighted average of partial derivative concerning spatial variable. From this, we obtain the system of nonlinear ordinary differential equations and it is linearized by using the quasilinearization method. Then by using the polynomial-based differential quadrature method for approximating derivative involving with spatial variable at specified grid point, we obtain the system of linear equation. Then they obtained linear system equation is solved by using the LU matrix decomposition method. To validate the applicability of the proposed method, two model examples are considered and solved at each specific grid point on its solution domain. The stability and convergent analysis of the present method is worked by supported the theoretical and mathematical statements and the accuracy of the solution is obtained. The accuracy of the present method has been shown in the sense of root mean square error norm and maximum absolute error norm and the local behavior of the solution is captured exactly. Numerical versus exact solutions and behavior of maximum absolute error between them have been presented in terms of graphs and the corresponding root mean square error norm and maximum absolute error norm presented in tables. The present method approximates the exact solution very well and it is quite efficient and practically well suited for solving the non-linear parabolic equation. The numerical result presented in tables and graphs indicates that the approximate solution is in good agreement with the exact solution.

Common Factors in Fraction-Free Matrix Decompositions

We consider LU and QR matrix decompositions using exact computations. We show that fraction-free Gauß–Bareiss reduction leads to triangular matrices having a non-trivial number of common row factors. We identify two types of common factors: systematic and statistical. Systematic factors depend on the reduction process, independent of the data, while statistical factors depend on the specific data. We relate the existence of row factors in the LU decomposition to factors appearing in the Smith–Jacobson normal form of the matrix. For statistical factors, we identify some of the mechanisms that create them and give estimates of the frequency of their occurrence. Similar observations apply to the common factors in a fraction-free QR decomposition. Our conclusions are tested experimentally.

Space charge effects and internal standardization in a four ion lenses interface: What is changing in quadrupole inductively coupled plasma mass spectrometry?

Space charge effects and internal standardization were evaluated in quadrupole inductively coupled plasma mass spectrometry (ICP-MS). Changes in analyte ion sensitivities as a function of analyte ion mass (11B+, 24Mg+, 45Sc+, 60Ni+, 71Ga+, 75As+, 111Cd+, 138Ba+ and 172Yb+) due to the presence of seven different matrix elements (23Na+, 63Cu+, 89Y+, 115In+, 159Tb+, 175Lu+ and 205Tl+) were measured by a current ICP-MS instrument arrangement. Elements with different masses (light, middle and high mass range) were evaluated as internal standard (IS): 7Li+, 88Sr+ and 209Bi+. Space charge effects generally becomes more severe for matrices containing elements with middle mass (115In+) and high masses (159Tb+, 175Lu+ and 205Tl+) but are not dependent on the masses of most analyte ions, differently from which was known about space charge effects when using first generation ICP-MS instruments. Thus, a single IS was effective for minimize the matrix induced changes (range from −20 to +20%) for analytes with different masses in In matrix (except to mass-to-charge ratio lower than 50); Tb matrix (except to 7Li+); Lu matrix (except to 175Yb+); and Tl matrix (except to 25Mg+).

Mixed Precision Block Fused Multiply-Add: Error Analysis and Application to GPU Tensor Cores

Computing units that carry out a fused multiply-add (FMA) operation with matrix arguments, referred to as tensor units by some vendors, have great potential for use in scientific computing. However, these units are inherently mixed precision, and existing rounding error analyses do not support them. We consider a mixed precision block FMA that generalizes both the usual scalar FMA and existing tensor units. We describe how to exploit such a block FMA in the numerical linear algebra kernels of matrix multiplication and LU factorization and give detailed rounding error analyses of both kernels. An important application is to GMRES-based iterative refinement with block FMAs, about which our analysis provides new insight. Our framework is applicable to the tensor core units in the NVIDIA Volta and Turing GPUs. For these we compare matrix multiplication and LU factorization with TC16 and TC32 forms of FMA, which differ in the precision used for the output of the tensor cores. Our experiments on an NVDIA V100 GPU confirm the predictions of the analysis that the TC32 variant is much more accurate than the TC16 one, and they show that the accuracy boost is obtained with almost no performance loss.

Real-Time HIL Emulation of Faulted Electric Machines Based on Nonlinear MEC Model

In electric machine drive systems, hardware-in-the-loop (HIL) emulation provides accurate testing of actual control system prototypes and protection devices interfaced with the electric machine model on a real-time simulator in a non-destructive environment particularly when faults are studied. A compromise between the model accuracy and computational burden makes the magnetic equivalent circuit (MEC) model ideal for real-time simulation of electric machines. However, satisfying the timing constraints of real-time simulation to accommodate internal machine faults is still challenging due to the nonlinearity and rotation of electric machines. In this paper, the transmission line modeling (TLM) method is utilized to keep the MEC coefficient matrix unchanged during nonlinear iterations. Afterward, for the first time, the entire potential of the TLM method for pre-calculation is exploited by proposing an efficient matrix re-ordering combined with the left-looking Gilbert–Peierls algorithm to minimize the computational burden of the sparse MEC matrix LU decomposition required in each time-step due to rotation. Furthermore, the massive hardware architecture of the field programmable gate array is used as the platform for implementation to fully exploit parallelism. With the proposed MEC-based real-time TLM method, the minimum time-step as low as 500 $\mu \text{s}$ can be achieved and the results validation with two-dimensional finite element model (FEM) of the commercial Jmag-Designer software shows the accuracy and efficiency of the proposed methodology.

An Algorithm to Calculate the Inverse Matrix if an Infinite Matrix

Abstract The purpose of this paper is to give an algorithm to calculate the left inverse matrix of an infinite matrix using the extension to infinite matrices of the LU matrix factorization. We give an example for the numerical solution of the heat equation.

Parallel real-world LU decomposition: Gauss vs. Crout algorithm

Abstract This paper presents numerical experiments with assorted versions of parallel LU matrix decomposition algorithms (Gauss and Crout algorithm). The tests have been carried out on the hardware platform with fourcore Skylake processor featuring hyperthreading technology doubling virtually core number. Parallelization algorithms have been implemented with the aid of classic POSIX threads library. Experiments have shown that basic 4-thread acceleration of all parallel implementations is almost equal to the number of threads/processors. Both algorithms are worth considering in real-world applications (Florida University collection). Gauss algorithm is a better performer, with respect to timing, in the case of matrices with lower density of nonzeros, as opposed to higher density matrices. The latter are processed more efficiently with the aid of Crout algorithm implementation.

An Efficient Load Flow Algorithm for AC–DC Distribution Systems

The AC–DC distribution systems have recently gained huge popularity due to advancements in power converters, high penetration of renewable energy resources and widespread use of DC loads. However, load flow in such systems is a challenging task, due to non-linear characteristics of power converters. This paper presents a new and computationally efficient load flow algorithm for AC distribution systems using graphical search technique along with matrix algebra and has been further extended for solving AC–DC distribution systems load flow problem by making use of proposed per-unit equivalent model of various power converter devices. In contrast with traditional methods, the proposed technique does not require any LU decomposition, matrix inversion and forward backward substitution of admittance matrix or Jacobian matrix. The significant feature of the proposed algorithm lies in formulation of path impedance matrix and loads beyond branch matrix, which will remain unaltered for the entire operation. Reconfiguration of these matrices takes place automatically in accordance with the change in input data. Because of the aforementioned reasons, the developed technique is computationally potent when applied on large sized distribution systems. The proposed method has been validated on various standard test systems. Results show feasibility and authenticity of the proposed method.

Minimum Breakpoint Set Determination for Directional Overcurrent Relay Coordination in Large-Scale Power Networks via Matrix Computations

Determination of the minimum breakpoint set (MBPS) is essential for optimum coordination of directional protection relays. However, in modern large-scale power networks, this task poses serious difficulties and becomes time consuming. This paper addresses such deficiencies and proposes a numerically efficient algorithm. Despite being derived from graph theory, the proposed algorithm functions via a series of matrix computations, that is, reduction, permutation, and LU factorization of the branch to bus incidence matrix of a given network. At first, parallel branches are united via the reduction step. Then, branches and buses are reset to form spanning trees of the network via the permutation step. Finally, MBPS is determined via the LU factorization step. Experimental results obtained by implementing the proposed algorithm in various test systems are presented to illustrate its computational performance.

Error estimates for the summation of real numbers with application to floating-point summation

Standard Wilkinson-type error estimates of floating-point algorithms involve a factor $$\gamma _k:=k\mathbf {u}/(1-k\mathbf {u})$$ for $$\mathbf {u}$$ denoting the relative rounding error unit of a floating-point number system. Recently, it was shown that, for many standard algorithms such as matrix multiplication, LU- or Cholesky decomposition, $$\gamma _k$$ can be replaced by $$k\mathbf {u}$$ , and the restriction on k can be removed. However, the arguments make heavy use of specific properties of both the underlying set of floating-point numbers and the corresponding arithmetic. In this paper, we derive error estimates for the summation of real numbers where each sum is afflicted with some perturbation. Recent results on floating-point summation follow as a corollary, in particular error estimates for rounding to nearest and for directed rounding. Our new estimates are sharp and unveil the necessary properties of floating-point schemes to allow for a priori estimates of summation with a factor omitting higher order terms.

The Unicorn Runtime: Efficient Distributed Shared Memory Programming for Hybrid CPU-GPU Clusters

Programming hybrid CPU-GPU clusters is hard. This paper addresses this difficulty and presents the design and runtime implementation of Unicorn —a parallel programming model for hybrid CPU-GPU clusters. In particular, this paper proves that efficient distributed shared memory style programing is possible and its simplicity can be retained across CPUs and GPUs in a cluster, minus the frustration of dealing with race conditions. Further, this can be done with a unified abstraction, avoiding much of the complication of dealing with hybrid architectures. This is achieved with the help of transactional semantics (on shared global address spaces), deferred bulk data synchronization, workload pipelining and various communication and computation scheduling optimizations. We describe the said abstraction, our computation and communication scheduling system and report its performance on a few benchmarks like Matrix Multiplication , LU Decomposition and 2D FFT . We find that parallelization of coarse-grained applications like matrix multiplication or 2D FFT using our system requires only about 30 lines of C code to set up the runtime. The rest of the application code is regular single CPU/GPU implementation. This indicates the ease of extending parallel code to a distributed environment. The execution is efficient as well. When multiplying two square matrices of size $65,536\times 65,536$ , Unicorn achieves a peak performance of 7.88 TFlop/s when run over a cluster of 14 nodes with each node equipped with two Tesla M2070 GPUs and two 6-core Intel Xeon 2.67 GHz CPUs, connected over a 32 Gbps Infiniband network. In this paper, we also demonstrate that the Unicorn programming model can be efficiently used to implement high level abstractions like MapReduce. We use such an extension to implement PageRank and report its performance. For a sample web of 500 million web pages, our implementation completes a page rank iteration in about 18 seconds (on average) on a 14-node cluster.

Factorization and Inversion of a Million Matrices using GPUs: Challenges and Countermeasures

This paper presents new algorithmic approaches and optimization techniques for LU factorization and matrix inversion of millions of very small matrices using GPUs. These problems appear in many scientific applications including astrophysics and generation of block-Jacobi preconditioners. We show that, for very small problem sizes, design and optimization of GPU kernels require a mindset different from the one usually used when designing LAPACK algorithms for GPUs. Techniques for optimal memory traffic, register blocking, and tunable concurrency are incorporated in our proposed design. We also take advantage of the small matrix sizes to eliminate the intermediate row interchanges in both the factorization and inversion kernels. The proposed GPU kernels achieve performance speedups vs. CUBLAS of up to 6× for the factorization, and 14× for the inversion, using double precision arithmetic on a Pascal P100 GPU.

Application of combined Newton–Raphson method to large load flow models

Load flow is the most used calculations in power system operation planning. Renewable resources have caused that system operator has to do power flow analysis and/or contingency analysis as fast as possible in order to predict next step in power system control. LU Decomposition of Jacobian matrix remains the most computationally expensive task during Newton–Raphson iterative method. Computational time appears to be critical issue when load flow calculation is performed on large power system load flow models. In this case, Jacobian matrix LU decomposition should not be performed in iterations in which convergence rate is not violated but performed in iterations in which convergence rate drop below specified level. In other words, modified Newton–Raphson Method which eliminates the repeated Jacobian matrix LU decomposition and generic Newton–Raphson method are combined depending on convergence rate. The paper presents application of proposed combined Newton–Raphson method which is based on convergence rate control. Comparison of combined, Shamanskii, generic, and modified Newton–Raphson methods is carried out taking into consideration computational time and number of iterations required to achieve convergence of load flow models of various dimensions.