Cholesky Factorization Algorithm Research Articles

Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes

In recent years, there has been widespread adoption of machine learning-based approaches to automate the solving of partial differential equations (PDEs). Among these approaches, Gaussian processes (GPs) and kernel methods have garnered considerable interest due to their flexibility, robust theoretical guarantees, and close ties to traditional methods. They can transform the solving of general nonlinear PDEs into solving quadratic optimization problems with nonlinear, PDE-induced constraints. However, the complexity bottleneck lies in computing with dense kernel matrices obtained from pointwise evaluations of the covariance kernel, and its partial derivatives, a result of the PDE constraint and for which fast algorithms are scarce. The primary goal of this paper is to provide a near-linear complexity algorithm for working with such kernel matrices. We present a sparse Cholesky factorization algorithm for these matrices based on the near-sparsity of the Cholesky factor under a novel ordering of pointwise and derivative measurements. The near-sparsity is rigorously justified by directly connecting the factor to GP regression and exponential decay of basis functions in numerical homogenization. We then employ the Vecchia approximation of GPs, which is optimal in the Kullback-Leibler divergence, to compute the approximate factor. This enables us to compute ϵ \epsilon -approximate inverse Cholesky factors of the kernel matrices with complexity O ( N log d ⁡ ( N / ϵ ) ) O(N\log ^d(N/\epsilon )) in space and O ( N log 2 d ⁡ ( N / ϵ ) ) O(N\log ^{2d}(N/\epsilon )) in time. We integrate sparse Cholesky factorizations into optimization algorithms to obtain fast solvers of the nonlinear PDE. We numerically illustrate our algorithm’s near-linear space/time complexity for a broad class of nonlinear PDEs such as the nonlinear elliptic, Burgers, and Monge-Ampère equations. In summary, we provide a fast, scalable, and accurate method for solving general PDEs with GPs and kernel methods.

Weighted tensor Golub–Kahan–Tikhonov-type methods applied to image processing using a t-product

This paper discusses weighted tensor Golub–Kahan-type bidiagonalization processes using the t-product. This product was introduced in M.E. Kilmer and C.D. Martin (2011). A few steps of a bidiagonalization process with a weighted least squares norm are carried out to reduce a large-scale linear discrete ill-posed problem to a problem of small size. The weights are determined by symmetric positive definite (SPD) tensors. Tikhonov regularization is applied to the reduced problem. An algorithm for tensor Cholesky factorization of SPD tensors is presented. The data is a laterally oriented matrix or a general third order tensor. The use of a weighted Frobenius norm in the fidelity term of Tikhonov minimization problems is appropriate when the noise in the data has a known covariance matrix that is not the identity. We use the discrepancy principle to determine both the regularization parameter in Tikhonov regularization and the number of bidiagonalization steps. Applications to image and video restoration are considered.

Exactly Solving Sparse Rational Linear Systems via Roundoff-Error-Free Cholesky Factorizations

Exactly solving sparse symmetric positive definite (SPD) linear systems is a key problem in mathematics, engineering, and computer science. This paper derives two new sparse roundoff-error-free (REF) Cholesky factorization algorithms which exactly solve sparse SPD linear systems $A {x} = {b}$, where $A \in \mathbb{Q}^{n x n}$ and ${x}, {b} \in {Q}^{n x p}$. The key properties of these factorizations are that (1) they exclusively use integer-arithmetic and (2) in the bit-complexity model, they solve the linear system $A {x} = {b}$ in time proportional to the cost of the integer-arithmetic operations. Namely, the overhead related to data structures and ancillary operations (those not strictly required to perform the factorization) is subsumed by the cost of the integer-arithmetic operations that are essential/intrinsic to the factorization. Notably, to date our algorithms are the only exact algorithm for solving SPD linear systems with this asymptotically efficient complexity bound. Computationally, we show that the novel factorizations are faster than both sparse rational-arithmetic LDL and sparse exact LU factorization. Altogether, the derived sparse REF Cholesky factorizations present a framework to solve any rational SPD linear system exactly and efficiently.

Sum of Kronecker products representation and its Cholesky factorization for spatial covariance matrices from large grids

The sum of Kronecker products (SKP) representation for spatial covariance matrices from gridded observations and a corresponding adaptive-cross-approximation-based framework for building the Kronecker factors are investigated. The time cost for constructing an n-dimensional covariance matrix is O(nk2) and the total memory footprint is O(nk), where k is the number of Kronecker factors. The memory footprint under the SKP representation is compared with that under the hierarchical representation and found to be one order of magnitude smaller. A Cholesky factorization algorithm under the SKP representation is proposed and shown to factorize a one-million dimensional covariance matrix in under 600 seconds on a standard scientific workstation. With the computed Cholesky factor, simulations of Gaussian random fields in one million dimensions can be achieved at a low cost for a wide range of spatial covariance functions.

Dynamic Scheduling Strategy for Block Parallel Cholesky Factorization Based on Activity on Edge Network

The efficient development of system software and design applications in parallel architecture is a notable challenge considering various aspects, such as load balancing, memory spaces, communication, and synchronization. This paper presents a block parallel Cholesky factorization algorithm for a multicore system, which is developed based on activity on edge network. First, the basic block computing tasks and their dependencies are taken as vertices and edges, respectively, and a directed acyclic graph corresponding to the specific block parallel Cholesky factorization is generated. Next, each edge of the directed acyclic graph is assigned to a weight equal to the processing time of the initial vertex of the edge, and the directed acyclic graph becomes an activity on edge network with only one starting and one ending vertex. Finally, a queuing algorithm is designed for the basic block computing tasks according to the edge activity on edge network, and a dynamic scheduling strategy is developed for block parallel Cholesky factorization. The results of the experiments concerning the parallel execution time of the algorithm in multicore systems with different configurations demonstrate that the proposed algorithm has notable advantages compared with the traditional static scheduling algorithm, and it exhibits satisfactory load balancing, parallelism, and scalability capacities.

A Schur complement approach to preconditioning sparse linear least-squares problems with some dense rows

The effectiveness of sparse matrix techniques for directly solving large-scale linear least-squares problems is severely limited if the system matrix A has one or more nearly dense rows. In this paper, we partition the rows of A into sparse rows and dense rows (As and Ad) and apply the Schur complement approach. A potential difficulty is that the reduced normal matrix AsTAs is often rank-deficient, even if A is of full rank. To overcome this, we propose explicitly removing null columns of As and then employing a regularization parameter and using the resulting Cholesky factors as a preconditioner for an iterative solver applied to the symmetric indefinite reduced augmented system. We consider complete factorizations as well as incomplete Cholesky factorizations of the shifted reduced normal matrix. Numerical experiments are performed on a range of large least-squares problems arising from practical applications. These demonstrate the effectiveness of the proposed approach when combined with either a sparse parallel direct solver or a robust incomplete Cholesky factorization algorithm.

A Multithreaded Algorithm for Sparse Cholesky Factorization on Hybrid Multicore Architectures

We present a multithreaded method for supernodal sparse Cholesky factorization on a hybrid multicore platform consisting of a multicore CPU and GPU. Our algorithm can utilize concurrency at different levels of the elimination tree by using multiple threads in both the CPU and the GPU. The elimination tree is a tree data structure describing the workflow of the factorization. Our experiments results on a platform consisting of an Intel multicore processor along with an Nvidia GPU indicate a significant improvement in performance and energy over single-threaded supernodal algorithm.

Accelerating sparse Cholesky factorization on GPUs

Sparse direct factorization is a critical component of scientific computingGPU acceleration of it is difficult due to algorithmic irregularity and PCIe perf.We present an algorithm for GPU acceleration which resolves these issues.The algorithm shows 2x speedup vs. best alternative CPU-only perf.Aspects of this Cholesky case should be useful in other algorithms (LDL, LU, ... ). Sparse factorization is a fundamental tool in scientific computing. As the major component of a sparse direct solver, it represents the dominant computational cost for many analyses. For factorizations which involve sufficient dense math, the substantial computational capability provided by GPUs (Graphics Processing Units) can help alleviate this cost. However, for many other cases, the prevalence of small/irregular dense math and the relatively slow communication between the host and device over the PCIe bus, make it challenging to significantly accelerate sparse factorization using the GPU.In this paper we describe a left-looking supernodal Cholesky factorization algorithm which permits improved utilization of the GPU when factoring sparse matrices. The central idea is to stream subtrees of the elimination tree through the GPU and perform the factorization of each subtree entirely on the GPU. This avoids the majority of the PCIe communication without the need for a complex task scheduler. Importantly, within these subtrees, many independent, small, dense operations are batched to minimize kernel launch overhead and many of these batched kernels are executed concurrently to maximize device utilization.Performance results for commonly studied matrices are presented along with suggested actions for further optimization.

Fail-Stop Failure Algorithm-Based Fault Tolerance for Cholesky Decomposition

Cholesky decomposition is a widely used algorithm to solve linear equations with symmetric and positive definite coefficient matrix. With large matrices, this often will be performed on high performance supercomputers with a large number of processors. Assuming a constant failure rate per processor, the probability of a failure occurring during the execution increases linearly with additional processors. Fault tolerant methods attempt to reduce the expected execution time by allowing recovery from failure. This paper presents an analysis and implementation of a fault tolerant Cholesky factorization algorithm that does not require checkpointing for recovery from fail-stop failures. Rather, this algorithm uses redundant data added in an additional set of processes. This differs from previous works with algorithmic methods as it addresses fail-stop failures rather than fail-continue cases. The proposed fault tolerance scheme is incorporated into ScaLAPACK and validated on the supercomputer Kraken. Experimental results demonstrate that this method has decreasing overhead in relation to overall runtime as the matrix size increases, and thus shows promise to reduce the expected runtime for Cholesky factorizations on very large matrices.

An efficient computational approach for evaluating radiation flux for laser driven inertial confinement fusion targets

Radiation flux computation on the target is very important for laser driven Inertial Confinement Fusion, and view-factor based equation models (MacFarlane, 2003; Srivastava et al., 2000) are often used to compute this radiation flux on the capsule or samples inside the hohlraum. However, the equation models do not lead to sparse matrices and may involve an intensive solution process when discrete mesh elements become smaller and the number of equations increases.An efficient approach for the computation of radiation flux is proposed in this paper, in which, (1) symmetric and positive definite properties are achieved by transformation, and (2) an efficient Cholesky factorization algorithm is applied to significantly accelerate such equations models solving process.Finally, two targets on a laser facility built in China are considered to validate the computing efficiency of present approach. The results show that the radiation flux computation can be accelerated by a factor of 2.

New rigorous perturbation bounds for the generalized Cholesky factorization

Some new rigorous perturbation bounds for the generalized Cholesky factorization with normwise or componentwise perturbations in the given matrix are obtained, where the componentwise perturbation has the form of backward rounding error for the generalized Cholesky factorization algorithm. These bounds can be much tighter than some existing ones while the conditions for them to hold are simple and moderate.

Convergence analysis of an algorithm for accurate inverse Cholesky factorization

This paper is concerned with factorization of symmetric and positive definite matrices which are extremely ill-conditioned. Following the results by Rump (1990), Oishi et al. (2007, 2009) and Ogita (2010), Ogita and Oishi (2012) derived an iterative algorithm for an accurate inverse matrix factorization based on Cholesky factorization for such ill-conditioned matrices. We analyze the behavior of the algorithm in detail and give reasons for convergency by the use of numerical error analysis. Main analysis is that each iteration reduces the condition number of a preconditioned matrix by a factor around the relative rounding error unit until convergence. This behavior is consistent with the numerical results.

Comparative Study on Zero Forcing Block Linear Equalization for TD-SCDMA

3G standards such as TD-SCDMA use multi-user detection (MUD) at the base station to combat multi address interference (MAI) and inter-symbol interference (ISI). We review Cholesky factorization and Schur algorithm for Zero Forcing Block Linear Equalization. We compare the solutions of different algorithm for ZF-BLE by computation, and analyze the complexity of those methods.

A modified algorithm for accurate inverse Cholesky factorization

A modified algorithm for accurate inverse Cholesky factorization

Solving a large scale radiosity problem on GPU-based parallel computers

The radiosity equation has been used widely in computer graphics and thermal engineering applications. The equation is simple to formulate but is challenging to solve when the number of Lambertian surfaces associated with an application becomes large. In this paper, we present the algorithms to compute the view factors and solve the set of radiosity equations using an out-of-core Cholesky decomposition method. This work details the algorithmic procedures of the computation of the view factors and the Cholesky solver. The data layout of the radiosity matrix follows the block cyclic decomposition scheme used in ScaLAPACK. The parallel computation of the view factors on the GPUs extends the algorithms based on a serial community code called view3d. To handle large matrices that exceed the device memory on GPU, an out-of-core algorithm for parallel Cholesky factorization is implemented. A performance study conducted on Keeneland, a hybrid CPU/GPU cluster at the National Institute for Computational Sciences, composed of 264 nodes of multicore CPU and GPU are shown and discussed.

A modified algorithm for accurate inverse Cholesky factorization

A modified algorithm for accurate inverse Cholesky factorization

Supernodal sparse Cholesky factorization on graphics processing units

SUMMARYSparse Cholesky factorization is the most computationally intensive component in solving large sparse linear systems and is the core algorithm of numerous scientific computing applications. A large number of sparse Cholesky factorization algorithms have previously emerged, exploiting architectural features for various computing platforms. The recent use of graphics processing units (GPUs) to accelerate structured parallel applications shows the potential to achieve significant acceleration relative to desktop performance. However, sparse Cholesky factorization has not been explored sufficiently because of the complexity involved in its efficient implementation and the concerns of low GPU utilization.In this paper, we present a new approach for sparse Cholesky factorization on GPUs. We present the organization of the sparse matrix supernode data structure for GPU and propose a queue‐based approach for the generation and scheduling of GPU tasks with dense linear algebraic operations. We also design a subtree‐based parallel method for multi‐GPU system. These approaches increase GPU utilization, thus resulting in substantial computational time reduction.Comparisons are made with the existing parallel solvers by using problems arising from practical applications. The experiment results show that the proposed approaches can substantially improve sparse Cholesky factorization performance on GPUs. Relative to a highly optimized parallel algorithm on a 12‐core node, we were able to obtain speedups in the range 1.59× to 2.31× by using one GPU and 1.80× to 3.21× by using two GPUs. Relative to a state‐of‐the‐art solver based on supernodal method for CPU‐GPU heterogeneous platform, we were able to obtain speedups in the range 1.52× to 2.30× by using one GPU and 2.15× to 2.76× by using two GPUs. Concurrency and Computation: Practice and Experience, 2013. Copyright © 2013 John Wiley & Sons, Ltd.

An Improved Square-Root Algorithm for V-BLAST Based on Efficient Inverse Cholesky Factorization

A fast algorithm for inverse Cholesky factorization is proposed, to compute a triangular square-root of the estimation error covariance matrix for Vertical Bell Laboratories Layered Space-Time arch...

Distributed SBP Cholesky factorization algorithms with near-optimal scheduling

The minimal block storage Distributed Square Block Packed (DSBP) format for distributed memory computing on symmetric and triangular matrices is presented. Three algorithm variants (Basic, Static, and Dynamic) of the blocked right-looking Cholesky factorization are designed for the DSBP format, implemented, and evaluated. On our target machine, all variants outperform standard full-storage implementations while saving almost half the storage. Communication overhead is shown to be virtually eliminated by the Static and Dynamic variants, both of which take advantage of hardware parallelism to hide communication costs. The Basic variant is shown to yield comparable or slightly better performance than the full-storage ScaLAPACK routine PDPOTRF while clearly outperformed by both Static and Dynamic. Models of execution assuming zero communication costs and overhead are developed. For medium- and larger-sized problems, the Static schedule is near optimal on our target machine based on comparisons with these models and measurements of synchronization overhead.

Out-of-Core Implementations of Cholesky Factorization: Loop-Based versus Recursive Algorithms

We compare, in the same framework, out-of-core implementations of the Cholesky factorization algorithm. The candidate implementations are the classical blocked left-looking variant and a more recent recursive formulation. Both have been implemented for real positive definite matrices: the former in the parallel out-of-core linear algebra package (POOCLAPACK) library and the latter in the scalable out-of-core linear algebra computations (SOLAR) library. We perform a theoretical analysis of the amount of input/output (I/O) operations required by each variant. We consider alternatives for the left-looking algorithm: the one-tile and two-tiles approaches. We show that when main memory is restricted, the one-tile approach yields less I/O volume. We then show that the left-looking implementation requires less I/O volume than the recursive variant. We have implemented all for complex matrices, and we report on numerical experiments.