In this paper, we propose some efficient data redistribution algorithms for redistributing matrices from 1D or 2D irregular format to block cyclic data distribution (BCDD) format, which can be much faster than the BLACS routinePXGEMR2D. These algorithms can be used to combine direct methods with iterative methods. The proposed algorithms divide the communication into two phases: one for processes in the same column and the other for processes in the same row, and the whole data redistribution task is divided into several independent sub-communications. The communication time can be reduced a lot compared with BLACS. Performance results show that our algorithms can be$2\times$–$5\times$faster than the BLACS routinePXGEMR2Dwhen using 4096 processes and the experiments are performed on Tianhe-2A supercomputer.

In this paper, we revisit the convergence of the Heavy-ball method, and present improved convergence complexity results in the convex setting. We provide the first non-ergodic O(1/k) rate result of the Heavy-ball algorithm with constant step size for coercive objective functions. For objective functions satisfying a relaxed strongly convex condition, the linear convergence is established under weaker assumptions on the step size and inertial parameter than made in the existing literature. We extend our results to multi-block version of the algorithm with both the cyclic and stochastic update rules. In addition, our results can also be extended to decentralized optimization, where the ergodic analysis is not applicable.

In this paper, we consider a class of nonconvex problems with linear constraints appearing frequently in the area of image processing. We solve this problem by the penalty method and propose the iteratively reweighted alternating minimization algorithm. To speed up the algorithm, we also apply the continuation strategy to the penalty parameter. A convergence result is proved for the algorithm. Compared with the nonconvex ADMM, the proposed algorithm enjoys both theoretical and computational advantages like weaker convergence requirements and faster speed. Numerical results demonstrate the efficiency of the proposed algorithm.

Nonconvex optimization algorithms with random initialization have attracted increasing attention recently. It has been showed that many first-order methods always avoid saddle points with random starting points. In this paper, we answer a question: can the nonconvex heavy-ball algorithms with random initialization avoid saddle points? The answer is yes! Direct using the existing proof technique for the heavy-ball algorithms is hard due to that each iteration of the heavy-ball algorithm consists of current and last points. It is impossible to formulate the algorithms as iteration like xk+1= g(xk) under some mapping g. To this end, we design a new mapping on a new space. With some transfers, the heavy-ball algorithm can be interpreted as iterations after this mapping. Theoretically, we prove that heavy-ball gradient descent enjoys larger stepsize than the gradient descent to escape saddle points to escape the saddle point. And the heavy-ball proximal point algorithm is also considered; we also proved that the algorithm can always escape the saddle point.

GEneral Matrix Multiply (GEMM) is the most fundamental computational kernel routine in the BLAS library. To achieve high performance, in-memory data must be prefetched into fast on-chip caches before they are used. Two techniques, software prefetching and data packing, have been used to effectively exploit the capability of on-chip least recent used (LRU) caches, which are popular in traditional high-performance processors used in high-end servers and supercomputers. However, the market has recently witnessed a new diversity in processor design, resulting in high-performance processors equipped with shared caches with non-LRU replacement policies. This poses a challenge to the development of high-performance GEMM in a multithreaded context. As several threads try to load data into a shared cache simultaneously, interthread cache conflicts will increase significantly. We present a Shared Cache Partitioning (SCP) method to eliminate interthread cache conflicts in the GEMM routines, by partitioning a shared cache into physically disjoint sets and assigning different sets to different threads. We have implemented SCP in the OpenBLAS library and evaluated it on Phytium 2000+, a 64-core AArch64 processor with private LRU L1 caches and shared pseudo-random L2 caches (per four-core cluster). Our evaluation shows that SCP has effectively reduced the conflict misses in both L1 and L2 caches in a highly optimized GEMM implementation, resulting in an improvement of its performance by 2.75% to 6.91%.

In image processing, total variation (TV) regularization models are commonly used to recover the blurred images. One of the most efficient and popular methods to solve the convex TV problem is the alternating direction method of multipliers (ADMM) algorithm, recently extended using the inertial proximal point method. Although all the classical studies focus on only a convex formulation, recent articles are paying increasing attention to the nonconvex methodology due to its good numerical performance and properties. In this paper, we propose to extend the classical formulation with a novel nonconvex alternating direction method of multipliers with the inertial technique (IADMM). Under certain assumptions on the parameters, we prove the convergence of the algorithm with the help of the Kurdyka-Łojasiewicz property. We also present numerical simulations on the classical TV image reconstruction problems to illustrate the efficiency of the new algorithm and its behavior compared with the well-established ADMM method.

Reconstructing the sparse signal from a few linear measurements has attracted increasing attentions in recent years. In this study, the authors propose the alternating projection (AP) method for sparse signal recovery with learning the sparsity of the original signal. Different with classical hard thresholding algorithms, the AP method regards the signal recovery problem as finding an intersect point of two sets. Theoretically, the authors prove that the proposed algorithm can reconstruct the s ‐sparse original signal provided the sensing matrix satisfies several assumptions when the noise is absent. They also prove that AP method is a noise‐robust algorithm, i.e. a tolerable reconstruction can be obtained by AP if the noise is small. In numerical experiments, the authors compare AP with several existing algorithms when being applied to sparse signals recovery and images reconstruction. The results demonstrate the efficiency of the proposed algorithm.

The nonsmooth and nonconvex regularization has many applications in imaging science and machine learning research due to its excellent recovery performance. A proximal iteratively reweighted nuclear norm algorithm has been proposed for the nonsmooth and nonconvex matrix minimizations. In this paper, we aim to investigate the convergence of the algorithm. With the Kurdyka-Łojasiewicz property, we prove the algorithm globally converges to a critical point of the objective function. The numerical results presented in this paper coincide with our theoretical findings.

This article describes a Matlab implementation of a mixed extended-precision package with three multiple-component formats. In the package, double-double, triple-double, and quad-double numbers, which are unevaluated sums of two, three and four IEEE double precision numbers respectively, are defined as Matlab classes including real and complex formats. We implement the four basic operations with mixed extended-precision based on these classes, and apply it to other operations and elementary functions. All operations and functions we presented are overloaded in Matlab. We experimentally evaluate this package that it allows the user to obtains more accurate numerical results with Matlab.

Polynomials are widely used in scientific computing and engineering. In this paper, we present an accurate and fast compensated algorithm to evaluate bivariate polynomials with floating-point coefficients. This algorithm is applying error free transformations to the bivariate Horner scheme and sum the final decomposition accurately. We also prove the proposed algorithm's accuracy with forward error analysis that the accuracy of the computed result is similar to the result computed by the bivariate Horner scheme in twice the working precision. Numerical experiments illustrate the behavior and it has higher efficiency than the bivariate Horner scheme implemented in double-double library.

Designing fast singular value decomposition (SVD) is significantly interesting in applications. The random direct SVD (RSVD) has provided a fast scheme to compute the well-approximate SVD by unilateral randomized sampling. In this paper, we present an efficient random algorithm in a bilateral sampling way. We also prove that the proposed algorithms can be bounded well and have less computational complexity compared to RSVD when the objective matrix is approximately square. Numerical experiments on graph Laplacian and Hilbert matrix demonstrate the efficiency and stability of the proposed methods.

This paper presents an implementation of an accurate and efficient compensated Double-precision General Matrix Multiplication (DGEMM) based on OpenBLAS for 64-bit ARMv8 multi-core processors. Due to cancellation phenomena in floating point arithmetic, the results of DGEMM may not be as accurate as expected. In order to increase the accuracy of DGEMM, we compensate the error introduced by its dot product kernel (GEBP) by applying an error-free transformation to rewrite the kernel in assembly language. We optimize the computations in the inner kernel through exploiting loop unrolling, instruction scheduling and software-implemented register rotation to exploit instruction level parallelism (ILP). We also conduct a priori error analysis of the derived CompDGEMM. Our compensated DGEMM is as accurate as the existing quadruple precision GEMM using MBLAS, but is up to 6.4x faster. Our parallel implementation achieves good performance and scalability under varying thread counts across a range of matrix sizes evaluated.

This paper presents the design and implementation of a highly efficient Double-precision General Matrix Multiplication (DGEMM) based on Open BLAS for 64-bit ARMv8 eight-core processors. We adopt a theory-guided approach by first developing a performance model for this architecture and then using it to guide our exploration. The key enabler for a highly efficient DGEMM is a highly-optimized inner kernel GEBP developed in assembly language. We have obtained GEBP by (1) maximizing its compute-to-memory access ratios across all levels of the memory hierarchy in the ARMv8 architecture with its performance-critical block sizes being determined analytically, and (2) optimizing its computations through exploiting loop unrolling, instruction scheduling and software-implemented register rotation and taking advantage of A64 instructions to support efficient FMA operations, data transfers and prefetching. We have compared our DGEMM implemented in Open BLAS with another implemented in ATLAS (also in terms of a highly-optimized GEBP in assembly). Our implementation outperforms the one in ALTAS by improving the peak performance (efficiency) of DGEMM from 3.88 Gflops (80.9%) to 4.19 Gflops (87.2%) on one core and from 30.4 Gflops (79.2%) to 32.7 Gflops (85.3%) on eight cores. These results translate into substantial performance (efficiency) improvements by 7.79% on one core and 7.70% on eight cores. In addition, the efficiency of our implementation on one core is very close to the theoretical upper bound 91.5% obtained from micro-benchmarking. Our parallel implementation achieves good performance and scalability under varying thread counts across a range of matrix sizes evaluated.

HPCG has become a new metric for the design and ranking of HPC. By incorporating a local symmetric Gauss-Seidel preconditioned, HPCG implements the Conjugate Gradient method to solve a sparse linear system. HPCG performs poorly with irregular memory access and may consume a great deal of MPI resources when it is executed on supercomputers. This paper focuses on optimizing SpMV and the Gauss-Seidel preconditioned, the two most important kernels in HPCG. By evaluating the performance impacts of several representative sparse matrix formats, ELLPACK is selected due to its suitability for SIMD, resulting in a speedup of 2.3x for the SpMV kernel. Multi-coloring is performed for Gauss-Seidel, resulting in a speedup of 7.3x over the reference implementation. The CG convergence rate may also be improved after multi-coloring. Our experimental results show that our optimization process works well on supercomputers, achieving 6.5 Gflops on a CPU-only node. This has boosted the total HPCG Gflops by about 7x, giving rise to 80,151 Gflops on 8192 CPU-only Tianhe-2 nodes.

This paper is concerned with the fast and accurate evaluation of elementary symmetric functions. We present a new compensated algorithm by applying error-free transformations to improve the accuracy of the so-called Summation Algorithm, which is used, by example, in the MATLAB's poly function. We derive a forward round off error bound and running error bound for our new algorithm. The round off error bound implies that the computed result is as accurate as if computed with twice the working precision and then rounded to the current working precision. The running error analysis provides a shaper bound along with the result, without increasing significantly the computational cost. Numerical experiments illustrate that our algorithm runs much faster than the algorithm using the classic double-double library while sharing similar error estimates. Such an algorithm can be widely applicable for example to compute characteristic polynomials from eigen values. It can also be used into the Rasch model in psychological measurement.