There are synergies of research interests and industrial efforts in modeling fairness and correcting algorithmic bias in machine learning. In this paper, we present a scalable algorithm for spectral clustering (SC) with group fairness constraints. Group fairness is also known as statistical parity where in each cluster, each protected group is represented with the same proportion as in the entirety. While FairSC algorithm (Kleindessner et al., 2019) is able to find the fairer clustering, it is compromised by high computational costs due to the algorithm deflation such that the resulting algorithm, called s-FairSC, only involves the sparse matrix-vector products and is able to fully exploit the sparsity of the fair SC model. The experimental results on the modified stochastic block model demonstrate that while it is comparable with FairSC in recovering fair clustering, s-FairSC is 12$\times$ faster than FairSC for moderate model sizes. s-FairSC is further demonstrated to be scalable in the sense that the computational costs of s-FairSC only increase marginally compared to the SC without fairness constraints.

We propose an efficient algorithm for solving orthogonal canonical correlation analysis (OCCA) in the form of trace-fractional structure and orthogonal linear projections. Even though orthogonality has been widely used and proved to be a useful criterion for visualization, pattern recognition and feature extraction, existing methods for solving OCCA problem are either numerically unstable by relying on a deflation scheme, or less efficient by directly using generic optimization methods. In this paper, we propose an alternating numerical scheme whose core is the sub-maximization problem in the trace-fractional form with an orthogonality constraint. A customized self-consistent-field (SCF) iteration for this sub-maximization problem is devised. It is proved that the SCF iteration is globally convergent to a KKT point and that the alternating numerical scheme always converges. We further formulate a new trace-fractional maximization problem for orthogonal multiset CCA and propose an efficient algorithm with an either Jacobi-style or Gauss-Seidel-style updating scheme based on the SCF iteration. Extensive experiments are conducted to evaluate the proposed algorithms against existing methods, including real-world applications of multi-label classification and multi-view feature extraction. Experimental results show that our methods not only perform competitively to or better than the existing methods but also are more efficient.

Some scientific and engineering applications need to compute a large number of eigenpairs of a large Hermitian matrix. Though the Lanczos method is effective for computing a few eigenvalues, it can be expensive for computing a large number of eigenpairs (e.g., in terms of computation and communication). To improve the performance of the method, in this paper, we study an s-step variant of thick-restart Lanczos (TRLan) combined with an explicit external deflation (EED). The s-step method generates a set of s basis vectors at a time and reduces the communication costs of generating the basis vectors. We then design a specialized matrix powers kernel (MPK) that reduces both the communication and computational costs by taking advantage of the special properties of the deflation matrix. We conducted numerical experiments of the new TRLan eigensolver using synthetic matrices and matrices from electronic structure calculations. The performance results on the Cori supercomputer at the National Energy Research Scientific Computing Center (NERSC) demonstrate the potential of the specialized MPK to significantly reduce the execution time of the TRLan eigensolver. The speedups of up to 3.1× and 5.3× were obtained in our sequential and parallel runs, respectively.

FAST-GE is a generalized spectral method for constrained clustering [Cucuringu et al., AISTATS 2016]. It incorporates the must-link and cannot-link constraints into two Laplacian matrices and then minimizes a Rayleigh quotient via solving a generalized eigenproblem, and is considered to be simple and scalable. However, there are two unsolved issues. Theoretically, since both Laplacian matrices are positive semi-definite and the corresponding pencil is singular, it is not proven whether the minimum of the Rayleigh quotient exists and is equivalent to an eigenproblem. Computationally, the locally optimal block preconditioned conjugate gradient (LOBPCG) method is not designed for solving the eigenproblem of a singular pencil. In fact, to the best of our knowledge, there is no existing eigensolver that is immediately applicable. In this paper, we provide solutions to these two critical issues. We prove a generalization of Courant-Fischer variational principle for the Laplacian singular pencil. We propose a regularization for the pencil so that LOBPCG is applicable. We demonstrate the robustness and efficiency of proposed solutions for constrained image segmentation. The proposed theoretical and computational solutions can be applied to eigenproblems of positive semi-definite pencils arising in other machine learning algorithms, such as generalized linear discriminant analysis in dimension reduction and multisurface classification via eigenvectors.

The Hubbard Hamiltonian provides a theoretical framework for describing electron interactions of quantum many-body systems in condensed matter physics. Determinant Quantum Monte Carlo (DQMC) simulations of the Hubbard Hamiltonian have contributed greatly to understanding important properties of materials. Physical measurements such as superconductivity and magnetic susceptibility are based on selected entries of a large set of Green's functions. The computations of Green's functions are equivalent to computing selected blocks of the inverses of large p-cyclic matrices. The performance of the state-of-art algorithm for computing Green's functions is around 100 Gflops on a 12-core Intel "Ivy Bridge" processor. In this paper, we describe a fast selected inversion (FSI) algorithm for computing selected entries of Green's functions and present a parallel implementation using hybrid MPI/OpenMP programming. The FSI algorithm rests on three ideas: (1) applying a block cyclic reduction for a structure-preserving reduction, (2) computing the inverse of the reduced block p-cyclic matrix by a structured orthogonal factorization, (3) using the block entries of the inverse of the reduced block p-cyclic matrix as seeds to rapidly form the selected inversion in parallel. Performance results of the new FSI algorithm on Edison, National Energy Research Scientific Computing Center (NERSC)'s Cray XC30 supercomputer, show an 80% improvement to 180 Gflops on the Intel "Ivy Bridge" processor. The parallel applications of the FSI algorithm for computing selected entries of multiple Green's functions reach to 20 -- 30 Tflops on 100 compute nodes with 2400 cores. The preliminary results show that the FSI algorithm speeds up a full DQMC simulation of the Hubbard Hamiltonian by a factor of five, reducing from three and a half hours down to only forty minutes on the l2-core processor.

A variation of Kublanovskaya's nonlinear QR method for solving banded nonlinear eigenvalue problems is presented in this article. The new method is iterative and specifically designed for problems too large to use dense linear algebra techniques. For the unstructurally banded nonlinear eigenvalue problem, a new data structure is used for storing the matrices to keep memory and computational costs low. In addition, an algorithm is presented for computing several nearby nonlinear eigenvalues to already-computed ones. Finally, numerical examples are given to show the efficacy of the new methods, and the source code has been made publicly available.

Numerical code uses floating-point arithmetic and necessarily suffers from roundoff and truncation errors. Error analysis is the process to quantify such uncertainty in the solution to a problem. Forward error analysis and backward error analysis are two popular paradigms of error analysis. Forward error analysis is more intuitive and has been explored and automated by the programming languages (PL) community. In contrast, although backward error analysis is more preferred by numerical analysts and the foundation for numerical stability, it is less known and unexplored by the PL community. To fill the gap, this paper presents an automated backward error analysis for numerical code to empower both numerical analysts and application developers. In addition, we use the computed backward error results to also compute the condition number, an important quantity recognized by numerical analysts for measuring how sensitive a function is to changes or errors in the input. Experimental results on Intel X87 FPU functions and widely-used GNU C Library functions demonstrate that our analysis is effective at analyzing the accuracy of floating-point programs.

The Determinant Quantum Monte Carlo (DQMC) method is one of the most powerful approaches for understanding properties of an important class of materials with strongly interacting electrons, including magnets and superconductors. It treats these interactions exactly, but the solution of a system of N electrons must be extrapolated to bulk values. Currently N 500 is state-of-the-art. Increasing N is required before DQMC can be used to model newly synthesized materials like functional multilayers. DQMC requires millions of linear algebra computations of order N matrices and scales as N3. DQMC cannot exploit parallel distributed memory computers efficiently due to limited scalability with the small matrix sizes and stringent procedures for numerical stability. Today, the combination of multisocket multicore processors and GPUs provides widely available platforms with new opportunities for DQMC parallelization. The kernel of DQMC, the calculation of the Green's function, involves long products of matrices. For numerical stability, these products must be computed using graded decompositions generated by the QR decomposition with column pivoting. The high communication overhead of pivoting limits parallel efficiency. In this paper, we propose a novel approach that exploits the progressive graded structure to reduce the communication costs of pivoting. We show that this method preserves the same numerical stability and achieves 70% performance of highly optimized DGEMM on a two-socket six-core Intel processor. We have integrated this new method and other parallelization techniques into QUEST, a modern DQMC simulation package. Using 36 hours on this Intel processor, we are able to compute accurately the magnetic properties and Fermi surface of a system of N = 1024 electrons. This simulation is almost an order of magnitude more difficult than N 500, owing to the N3 scaling. This increase in system size will allow, for the first time, the computation of the magnetic and transport properties of layered materials with DQMC. In addition, we show preliminary results which further accelerate DQMC simulations by using GPU processors.

A matrix operation is referred to as a hard-to-parallel matrix operation (HPMO) if it has serial bottlenecks that are hardly parallelizable. Otherwise, it is referred to as an easy-to-parallel matrix operation (EPMO). Empirical evidences showed the performance scalability of an HPMO is significantly poorer than an EPMO on multicore and distributed architectures. As the result, the design of higher-level algorithms for applications, for the performance considerations on multicore and distributed architectures, should avoid the use of HPMOs as the computational kernels. In this paper, as a case study, we present an HPMO-avoiding algorithm for the Green's function calculation in quantum Monte Carlo simulation. The original algorithm utilizes the QR-decomposition with column pivoting (QRP) as its computational kernel. QRP is an HPMO. The redesigned algorithm maintains the same simulation stability but employs the standard QR decomposition without pivoting (QR), which is an EPMO. Different implementations of the redesigned algorithm on multicore and distributed architectures are investigated. Although some implementations of the redesigned method use about a factor of three more floating-point operations than the original algorithm, they are about 20% faster on a quad core system and 2.5 times faster on a 1024-CPU massively parallel processing system. The broader impact of the redesign of higher-level matrix algorithms to avoid HPMOs in other computational science applications is also discussed.

Determinant Quantum Monte Carlo (DQMC) simulation has been widely used to reveal macroscopic properties of strong correlated materials. However, parallelization of the DQMC simulation is extremely challenging duo to the serial nature of underlying Markov chain and numerical stability issues. We extend previous work with novelty by presenting a hybrid granularity parallelization (HGP) scheme that combines algorithmic and implementation techniques to speed up the DQMC simulation. From coarse-grained parallel Markov chain and task decompositions to fine-grained parallelization methods for matrix computations and Green's function calculations, the HGP scheme explores the parallelism on different levels and maps the underlying algorithms onto different computational components that are suitable for modern high performance heterogeneous computer systems. Practical techniques, such as communication and computation overlapping, message compression and load balancing are also considered in the proposed HGP scheme. We have implemented the DQMC simulation with the HGP scheme on an IBM Blue Gene/P system. The effectiveness of the new scheme is demonstrated through both theoretical analysis and performance results. Experiments have shown over a factor of 80 speedups on an IBM Blue Gene/P system with 1,014 computational processors.

The Thick-Restart Lanczos (TRLan) method is an effective method for solving large-scale Hermitian eigenvalue problems. The performance of the method strongly depends on the dimension of the projection subspace used at each restart. In this article, we propose an objective function to quantify the effectiveness of the selection of subspace dimension, and then introduce an adaptive scheme to dynamically select the dimension to optimize the performance. We have developed an open-source software package a --TRLan to include this adaptive scheme in the TRLan method. When applied to calculate the electronic structure of quantum dots, a --TRLan runs up to 2.3x faster than a state-of-the-art preconditioned conjugate gradient eigensolver.

This paper presents a multiparameter moment-matching-based model order reduction technique for parameterized interconnect networks via a novel two-directional Arnoldi process (TAP). It is referred to as a Parameterized Interconnect Macromodeling via a TAP (PIMTAP) algorithm. PIMTAP inherits the advantages of previous multiparameter moment-matching algorithms and avoids their shortfalls. It is numerically stable and adaptive. PIMTAP model yields the same form of the original state equations and preserves the passivity of parameterized RLC networks like the well-known method passive reduced-order interconnect macromodeling algorithm for nonparameterized RLC networks.

We describe an efficient implementation and present a performance study of an automated multi-level substructuring (AMLS) method for sparse eigenvalue problems. We assess the time and memory requirements associated with the key steps of the algorithm, and compare it with the shift-and-invert Lanczos algorithm. Our eigenvalue problems come from two very different application areas: accelerator cavity design and normal-mode vibrational analysis of polyethylene particles. We show that the AMLS method, when implemented carefully, outperforms the traditional method in broad application areas when large numbers of eigenvalues are sought, with relatively low accuracy.

This paper presents a multiparameter moment- matching based model order reduction technique for parameter- ized interconnect networks via a novel two-directional Arnoldi process. It is referred to as a PIMTAP algorithm, which stands for Parameterized Interconnect Macromodeling algorithm via a Two-directional Arnoldi Process. PIMTAP inherits the advan- tages of previous multiparameter moment-matching algorithms and avoids their shortfalls. It is numerically stable and adaptive, and preserves the passivity of parameterized RLC networks.

Extracting features from point sets is becoming increasingly important for purposes like model classification, matching, and exploration. We introduce a technique for segmenting a point-sampled surface into distinct features without explicit construction of a mesh or other surface representation. Our approach achieves computational efficiency through a three-phase segmentation process. The first phase of the process uses a topological approach to define features and coarsens the input, resulting in a set of supernodes, each one representing a collection of input points. A graph cut is employed in the second phase to bisect the set of supernodes. Similarity between supernodes is computed as a weighted combination of geodesic distances and connectivity. Repeated application of the graph cut results in a hierarchical segmentation of the point input. In the last phase, a segmentation of the original point set is constructed by refining the segmentation of the supernodes based on their associated feature sizes.We apply our segmentation algorithm on laser-scanned models to evaluate its ability to capture geometric features in complex data sets.

Recently model order reduction techniques for second-order systems have obtained many research interests for the simulation of RCS interconnect circuits employing susceptance elements. In this paper, we propose a Block SAPOR (Block Second-order Arnoldi method for Passive Order Reduction) for Multi-Input Multi-Output RCS Circuits. The proposed Block SAPOR algorithm can simultaneously guarantee passivity and achieve higher accuracy than the first order reduction technique PRIMA. Most importantly, the reduced system matrices obtained by the proposed method can preserve the structure of the original system matrices. Such a nice property makes it possible to construct an equivalent RCS circuit for the reduced system.

The recently-introduced susceptance element exhibits many prominent features in modeling the on-chip magnetic couplings. For an RCS circuit, it is better to be formulated as a second-order system. Therefore, corresponding MOR (model-order reduction) techniques for second-order systems are desired to efficiently deal with the ever-increasing circuit scale and to preserve essential model properties. We first review the existing MOR methods for RCS circuits, such as ENOR and SMOR, and discuss several key issues related to numerical stability and accuracy of the methods. Then, a technique, SAPOR (second-order Arnoldi method for passive order reduction), is proposed to effectively address these issues. Based on an implementation of a generalized second-order Arnoldi method, SAPOR is numerically stable and efficient. Meanwhile, the reduced-order system also guarantees passivity.

Recently, there has been a great deal of interest in using the Pade Via Lanczos (PVL) technique to analyze the transfer functions and impulse responses of large-scale linear circuits. In this paper, a matrix-based derivation of the error between the original circuit transfer function and the reduced-order transfer function generated using the PVL technique is presented. This error measure may be used for the development of an automated termination of the Lanczos process in the PVL technique and achieve the desired accuracy of the approximate transfer function. PVL coupled with such an error bound will be referred to as the PVL-WEB algorithm.

SRRT is a Fortran program to calculate an approximate orthonomral basis fr a dominant invariant subspace of a real matrix A by the method of simultaneous iteration. Specifically, given an integer m , SRRIT computes a matrix Q with m orthonormal columns and real quasi-triangular matrix T or order m such that the equation AQ = QT is satisfied up to a tolerance specified by the user. The eigenvalues of T are approximations to the m eigenvalues of largest absolute magnitude of A and the columns of Q span the invariant subspace corresponding to those eigenvalues. SRRIT references A only through a user-provided subroutine to form the product AQ ; hence it is suitable for large sparse problems.

We review the theory of condition numbers for the nonsymmetric eigenproblem and give a tabular summary of bounds for eigenvalues, means of clusters of eigenvalues, eigenvectors, invariant subspaces, and related quantities. We describe the design of new algorithms for estimating these condition numbers. Fortran subroutines implementing these algorithms are in the LAPACK library [1].

The goal of the LAPACK project is to design and implement a portable linear algebra library for efficient use on a variety of high-performance computers. The library is based on the widely used LINPACK and EISPACK packages for solving linear equations, eigenvalue problems, and linear least-squares problems, but extends their functionality in a number of ways. The major methodology for making the algorithms run faster is to restructure them to perform block matrix operations (e.g., matrix-matrix multiplication) in their inner loops. These block operations may be optimized to exploit the memory hierarchy of a specific architecture. The LAPACK project is also working on new algorithms that yield higher relative accuracy for a variety of linear algebra problems.>