Search

With the advent of massive datasets, distributed techniques for processing information and carrying out computations are expected to enable exceptional possibilities for engineering, data intensive sciences, and better decision making. Unfortunately, distributive computation networks are prone to straggling servers. In recent years, utilizing coding-theoretic techniques has proven to be a very powerful tool for recovering either exact or approximate computations in the presence of stragglers. Furthermore, existing algorithms for mathematical programming, which are the core component of techniques for processing large datasets, often prove ineffective for scaling to the extent of all available data. Over the past two decades, randomized dimensionality reduction has proven to be a promising tool for approximate computations over large datasets. These two approaches to dealing with large datasets are referred to as ``Coded Computing’’, and ``Randomized Numerical Linear Algebra’’ or ``Sketching’’. In this dissertation, we introduce new methods for both these areas, extend existing methods, combine techniques from the two areas, and show how sketching algorithms can be utilized to devise coded computing schemes. Additionally, in certain schemes, we provide security through sketching and polynomial codes. The applications we focus on are distributed steepest descent and stochastic steepest descent, matrix-matrix multiplication, graph sparsification, linear regression, matrix inversion and pseudo-inversion, which are vital components to optimization, machine learning, and arguably all engineering disciplines.

Methodologies that ensure the compressibility of tensors are introduced. Bounds on the storage costs with respect to various tensor formats are derived. A new algorithm combining data-sparse tensor formats and factored alternating direction implicit method is designed to solve Sylvester tensor equations, and incorporated in a fast spectral Poisson equation solver on cubes with optimal complexity. New parallelizable algorithms for computing the tensor-train decomposition of tensors in original format, streaming data, Tucker format, and that satisfy algebraic relations, are proposed. Based on the input format, the algorithms involve deterministic or probabilistic aspects, and all have guarantees of accuracy. Scaling analysis and numerical experiments are provided to demonstrate computational and storage efficiency. An ultraspherical spectral method is developed for fractional partial differential equations via the Caffarelli--Silvestre extension on disk and rectangular domains. A parallel domain decomposition solver is designed for multi-core performance of non-smooth functions. The discretized equation is solved via direct tensor equation solvers, and numerical performance is shown with a fractional PDE constrained optimization problem. Linear systems in electron correlation calculation from computational chemistry are converted into tensor equations to reduce computing and storage costs. Several algorithms are developed to exploit the sparsity and data-sparsity of chemical structures. Numerical results indicate that tensor equation solvers are competitive over traditional linear system solvers with both canonical and localized orbital bases formulations. The quantized tensor-train format of tensors is introduced to approximate analytic functions via Chebyshev polynomial expansions. Analysis of different types of singularities is carried out, leading to theoretical guarantees of coefficient storage compressibility.

In this thesis we study the computational problem of approximately diagonalizing an arbitrary complex, n⨉n matrix A in floating point arithmetic — that is, finding V invertible and D diagonal with ||A - VDV⁻¹|| ≤ δ||A|| for some desired accuracy δ>0. Our major contributions are the following.(i) We prove that a randomized variant of the Spectral Bisection algorithm can approximately diagonalize any n⨉n matrix using O(Tₘₘ(n) log(n/δ) log n) arithmetic operations, on a floating point machine with O(log⁴(n/δ) log n) bits of precision, with probability 1 - O(1/n) — where Tₘₘ(n) is the arithmetic cost of numerically stable n⨉n matrix multiplication.(ii) We prove that for every k = 2,4,8,... and B ≥ 1 there is a shifting strategy for the Shifted QR algorithm, S(k,B), which converges rapidly in exact arithmetic on every Hessenberg matrix H₀ with eigenvector condition number less than B (e.g. of the form H₀ = VDV⁻¹ with D diagonal and ||V||||V⁻¹||≤B). The resulting algorithm gives a sequence H₀,H₁,H₂,... of Hesssenberg matrices unitarily similar to H₀, and rapid convergence means: each iteration costs O(n² k log k + B^(16k⁻¹ log k)) arithmetic operations, and for any ω>0, it takes O(log(1/ω)) iterations to drive the smallest subdiagonal entry below ω||H₀||.(iii) We prove that the shifting strategy S(k,B) can be implemented in floating point as a randomized algorithm, using O(k log(nB||H₀||/ωgap(H₀))) bits of precision (where gap(H) is the minimum eigenvalue gap, or the smallest distance between two eigenvalues of H), and succeeding with probability 1 - O(1/n). For any n⨉n matrix A, this can be used to give a randomized algorithm for computing the matrix D in the notion of approximate diagonalization above (e.g., the approximate eigenvalues a matrix δ||A||-close to A), using O(n³ log²(n/δ) poly(log(log(n/δ)))) arithmetic operations, on a floating point machine with O(log²(n/δ) log(log(n/δ))) bits of precision, again succeeding with probability 1-O(1/n).A crucial step along the way is a new result in random matrix theory: for any complex n⨉n matrix A and δ > 0, if G is an n⨉n matrix with independent, standard complex Gaussian entries, then with probability at least 1 - O(1/n), A + δG has eigenvector condition number at most n²/δ and minimum eigenvalue gap at least δ/n². We show similar results for random real perturbations of real matrices as well.

Recent years have seen rapid progress in the field of quantum algorithms for linear algebra problems. These algorithms can solve important problems in quantum chemistry, condensed matter physics, and quantum field theory simulation, with potentially exponential speedup compared to classical algorithms. The progress is in part due to the development of new methods to implement matrix functions on a quantum computer. Examples include the linear combination of unitaries (LCU) method, quantum signal processing (QSP), and quantum singular value transformation (QSVT), which is closely related to QSP. Using these methods, one can construct matrix functions to filter out unwanted eigenstates of a given Hermitian matrix, and thereby obtain a target eigenstate on a quantum computer. This technique we call quantum eigenstate filtering.We focus our attention on two specific problems: ground state preparation and energy estimation, and solving linear systems. For the ground state preparation problem, we want to prepare the ground state, i.e., the lowest eigenstate of the Hamiltonian $H$, while for ground state energy estimation, we want to estimate the ground state energy, i.e., the the lowest eigenvalue of the Hamiltonian $H$. The ground state energy estimation problem is cannot be solved in polynomial time even on a quantum computer in the worst case. However, this problem becomes tractable with additional assumptions, which may be satisfied in many real-world applications. In this work we assume access to a good initial guess of the ground state, such that its overlap with the ground state is lower bounded by a parameter $\gamma$. Under this assumption, we present quantum algorithms based on the quantum eigenstate filtering technique to estimate the ground state energy and to prepare the ground state, with nearly optimal scaling with respect to both $\gamma$ and the precision.The above mentioned algorithms work in the setting where the quantum computers are fully fault-tolerant, i.e., the noise is sufficiently suppressed through quantum error correction. We also present an algorithm to estimate the ground state energy in the early fault-tolerant setting, in which the the algorithm is limited by the number of qubits available and the coherence time. In this setting our algorithm uses a simple circuit and lower circuit depth at the expense of longer runtime. However the algorithm can still achieve the optimal Heisenberg-limited precision scaling despite these restrictions.The quantum eigenstate filtering technique is also useful in solving the quantum linear system problem, in which we are asked to solve $Ax=b$ on a quantum computer. We show that the problem can be solved with a query complexity almost matching its lower bound by traversing the adiabatic path using quantum eigenstate filtering.Chapter 1 introduces the basic problem setup and provides an overview of the results in this paper. Chapter 2 presents the nearly optimal algorithms for ground state preparation and energy estimation. Chapter 3 discussed implementation on early fault-tolerant devices. Chapter 4 applies quantum eigenstate filtering to solve the quantum linear system problem. Chapters 2, 3, and 4 are based on [Lin, Tong, Quantum 4 (2020), p.~372], [Lin, Tong, PRX Quantum 3.1 (2022), p.~010318], [Lin, Tong, Quantum 4 (2020), p.~361] respectively, all of which are published under CC BY 4.0.

New numerical methods using rational functions are presented for applications in linear algebra and signal processing. Classical results from Zolotarev are applied to develop a collection of low rank methods and theoretical results for computing with matrices that have special displacement structures using the alternating direction implicit (ADI) method. This includes a new low rank method for solving Sylvester and Lyapunov matrix equations with right hand sides that have decaying singular values, spectrally accurate low rank solvers for certain elliptic partial differential equations with smooth right-hand sides, and explicit bounds on the singular values of special families of structured matrices. Methods from conformal mapping and adaptive rational approximation are applied to build approximate solutions to Zolotarev's problem on sets where solutions are not known. This leads to new bounds on the numerical ranks of matrices, and it generalizes the regime in which ADI-based methods can be applied. The approximate solutions supply quasi-optimal ADI shift parameters for solving Sylvester matrix equations. A superfast rank-structured solver for Toeplitz linear systems is designed with ADI-based compression methods, and theoretical arguments are supplied that justify the effectiveness of rank-structured solvers for Toeplitz and related linear systems. The solvers are competitive with the state of the art, and rational approximation arguments are used to derive explicit error bounds on the numerical ranks of important submatrices for various weakly admissible hierarchical formats. A data-driven rational approximation framework is developed for reconstructing signals from samples with poorly separated spectral content. This approach combines a variant of Prony's method with a modified version of the AAA algorithm to construct representations of signals in both frequency and time space. The approximation methods are automatic and adaptive, requiring no tuning or manual parameter selection, and they are robust to various forms of corruption, including additive Gaussian noise, perturbed sampling grids, and missing data. A collection of algorithms and an accompanying software package for adaptively computing with these representations is introduced that includes procedures for differentiation/integration, rootfinding/polefinding, convolution, filtering, extrapolation, and more.

In this thesis, we present novel numerical algorithms in machine learning and data analysis. This thesis consists of two chapters. Each chapter is self-contained.In chapter 1, we present novel methods for performing k-fold cross validation for ridge regression. Ridge regression is a very broadly used method for reduce overfitting in linear regression. And it is used in many aspects of data analysis. One very important question for ridge regression is to find the best regularization hyperparameter $\lambda$. However, this can be a time consuming procedure. Here we present an efficient yet numerical stable algorithm for computing the relative error of ridge regression across different $\lambda$'s. Besides that, we provide an novel algorithm for finding the best hyper parameter $\lambda$ without computing all the relative errors.In chapter 2, we present algorithms for efficiently finding low rank approximation of a large sparse matrix with improved Lanczos Algorithm. Low rank approximation is one of the most important techniques in data analysis. In particular, people do principal component analysis based on that. We tackle this problem by finding the truncated SVD of the matrix using our improved version of Lanczos matrix. Our method guarantees that we only perform restart of the algorithm when necessary. Also, we prove a novel scheme for doing reorthogonalization during the Lanczos iteration. Besides that, we provide new stopping criteria for Lanczos algorithm directly based on the quality of the low rank approximation.

Mixed-precision algorithms are a class of algorithms that uses low precision in part of the algorithm in order to save time and energy with less accurate computation and communication. These algorithms usually utilize iterative refinement processes to improve the approximate solution obtained from low precision to the accuracy we desire from doing all the computation in high precision. Due to the demand of deep learning applications, there are hardware developments offering different low-precision formats including half precision (FP16), Bfloat16 and integer operations for quantized integers, which uses integers with a shared scalar to represent a set of equally spaced numbers. As new hardware architectures focus on bringing performance in these formats, the mixed-precision algorithms have more potential leverage on them and outmatch traditional fixed-precision algorithms. This dissertation consists of two articles. In the first article, we adapt one of the most fundamental algorithms in numerical linear algebra---LU factorization with partial pivoting--- to use integer arithmetic. With the goal of obtaining a low accuracy factorization as the preconditioner of generalized minimal residual (GMRES) to solve systems of linear equations, the LU factorization is adapted to use two different fixed-point formats for matrices L and U. A left-looking variant is also proposed for matrices with unbounded column growth. Finally, GMRES iterative refinement has shown that it can work on matrices with condition numbers up to 10000 with the algorithm that uses int16 as input and int32 accumulator for the update step. The second article targets symmetric and Hermitian eigenvalue problems. In this section we revisit the SICE algorithm from Dongarra et al. By applying the Sherman-Morrison formula on the diagonally-shifted tridiagonal systems, we propose an updated SICE-SM algorithm. By incorporating the latest two-stage algorithms from the PLASMA and MAGMA software libraries for numerical linear algebra, we achieved up to 3.6x speedup using the mixed-precision eigensolver with the blocked SICE-SM algorithm for iterative refinement when compared with full double complex precision solvers for the cases with a portion of eigenvalues and eigenvectors requested.

This dissertation is about computational tools based on randomized numerical linear algebra for handling larg-scale matrix data. Since large datasets have become commonly available in a wide variety of modern applications, there has been an increasing demand for numerical methods for storing, processing, and learning from them. Matrices, the classical form for representing datasets, naturally connect these tasks with the rich literature of numerical linear algebra. For a diverse collection of problems, randomized methods offer extraordinary efficiency and flexibility. This work focuses on using randomized numerical linear algebra to build practical algorithms for problems of massive size and high complexity that traditional methods are unable to handle. Through this dissertation, we explore topics across network science, Gaussian process regression, natural language processing, and quantum chemistry. Our contribution includes a collection of scalable and robust numerical methods under a unifying theme, accompanied by efficient implementations. As a result, we are able to significantly speed up the computation for several existing applications, and explore problems and datasets that were intractable before.

Probing methods for trace estimation of large, sparse matrices has been studied for several decades. In recent years, there has been some work to extend these techniques to instead estimate the diagonal entries of these systems directly. We extend some analysis of trace estimators to their corresponding diagonal estimators, propose a new class of deterministic diagonal estimators which are well-suited to parallel architectures along with heuristic arguments for the design choices in their construction, and conclude with numerical results on diagonal estimation and ordering problems, demonstrating the strengths of our newly-developed methods alongside existing methods.

In this dissertation, we study two problems that have the theme of extracting information from lower dimensional samples. A number of recent works have studied algorithms for entrywise $\ell_p$-\emph{low rank approximation}, namely algorithms which given an $n \times d$ matrix $A$ (with $n \geq d$), output a rank-$k$ matrix $B$ minimizing $\|A-B\|_p^p=\sum_{i,j} |A_{i,j} - B_{i,j}|^p$ when $p > 0$; and $\|A-B\|_0 = \sum_{i,j} [A_{i,j} \neq B_{i,j}]$ for $p=0$, where $\|A-B\|_0$ denotes the number of entries $(i,j)$ for which $A_{i,j} \neq B_{i,j}$.For $p = 1$, this is often considered more robust than the SVD, while for $p = 0$ this corresponds to minimizing the number of disagreements, or robust PCA. This problem is known to be NP-hard for $p \in \{0,1\}$, already for $k = 1$, and while there are polynomial time approximation algorithms, their approximation factor is at best $\poly(k)$. It was left open if there was a polynomial-time approximation scheme (PTAS) for $\ell_p$-approximation for any $p \geq 0$. We show the following:\begin{enumerate} \item On the algorithmic side, for $p \in (0,2)$, we use a technique called \emph{sketching} to give the first $n^{\poly(k/\eps)}$ time $(1+\eps)$-approximation algorithm. For $p = 0$, there are various problem formulations, a common one being the binary setting for which $A\in\{0,1\}^{n\times d}$ and $B = U \cdot V$, where $U\in\{0,1\}^{n \times k}$ and $V\in\{0,1\}^{k \times d}$. For this setting, we obtain an algorithm with time $n \cdot d^{\poly(k/\eps)}$. \item On the hardness front, for $p \in (1,2)$, we show under the Small Set Expansion Hypothesis and Exponential Time Hypothesis (ETH), there is no constant factor approximation algorithm running in time $2^{k^{\delta}}$ for a constant $\delta > 0$, showing an exponential dependence on $k$ is necessary. We also show for finite fields of constant size, under the ETH, that any fixed constant factor approximation algorithm requires $2^{k^{\delta}}$ time for a constant $\delta > 0$.\end{enumerate} \emph{Population recovery} is the problem of learning an unknown distribution over an unknown set of $n$-bit strings, given access to independent \emph{traces} from the distribution that have been independently corrupted according to some noise channel. Recent work has intensively studied such problems both for the {bit-flip} noise channel and for the erasure noise channel. In this dissertation we initiate the study of population recovery under the \emph{deletion channel}, in which each bit $b$ is independently \emph{deleted} with some fixed probability and the surviving bits are concatenated and transmitted. This is a far more challenging noise model than bit-flip~noise or erasure noise; indeed, even the simplest case in which the population is of size 1 (corresponding to a trivial probability distribution supported on a single string) corresponds to the \emph{trace reconstruction} problem, which is a challenging problem that has received much recent attention. In this work we give algorithms and lower bounds for population recovery under the deletion channel when the population size is some value $\ell > 1$. As our main sample complexity upper bound, we show that for any population size $\ell = o(\log n / \log \log n)$, a population of $\ell$ strings from $\zo^n$ can be learned under deletion channel noise using $\smash{2^{n^{1/2 + o(1)}}}$ samples. On the lower bound side, we show that at least $n^{\Omega{(\ell)}}$ samples are required to perform population recovery under the deletion channel when the population size is $\ell$, for all $\ell \leq n^{1/2 - \epsilon}$.

This thesis is devoted to answering a range of questions in machine learning and data analysis related to large dimensional matrices and scientific data. Two key research objectives connect the different parts of the thesis: (a) development of fast, efficient, and scalable algorithms for machine learning which handle large matrices and high dimensional data; and (b) design of learning algorithms for scientific data applications. The work combines ideas from multiple, often non-traditional, fields leading to new algorithms, new theory, and new insights in different applications. The first of the three parts of this thesis explores numerical linear algebra tools to develop efficient algorithms for machine learning with reduced computation cost and improved scalability. Here, we first develop inexpensive algorithms combining various ideas from linear algebra and approximation theory for matrix spectrum related problems such as numerical rank estimation, matrix function trace estimation including log-determinants, Schatten norms, and other spectral sums. We also propose a new method which simultaneously estimates the dimension of the dominant subspace of covariance matrices and obtains an approximation to the subspace. Next, we consider matrix approximation problems such as low rank approximation, column subset selection, and graph sparsification. We present a new approach based on multilevel coarsening to compute these approximations for large sparse matrices and graphs. Lastly, on the linear algebra front, we devise a novel algorithm based on rank shrinkage for the dictionary learning problem, learning a small set of dictionary columns which best represent the given data. The second part of this thesis focuses on exploring novel non-traditional applications of information theory and codes, particularly in solving problems related to machine learning and high dimensional data analysis. Here, we first propose new matrix sketching methods using codes for obtaining low rank approximations of matrices and solving least squares regression problems. Next, we demonstrate that codewords from certain coding scheme perform exceptionally well for the group testing problem. Lastly, we present a novel machine learning application for coding theory, that of solving large scale multilabel classification problems. We propose a new algorithm for multilabel classification which is based on group testing and codes. The algorithm has a simple inexpensive prediction method, and the error correction capabilities of codes are exploited for the first time to correct prediction errors. The third part of the thesis focuses on devising robust and stable learning algorithms, which yield results that are interpretable from specific scientific application viewpoint. We present Union of Intersections (UoI), a flexible, modular, and scalable framework for statistical-machine learning problems. We then adapt this framework to develop new algorithms for matrix decomposition problems such as nonnegative matrix factorization (NMF) and CUR decomposition. We apply these new methods to data from Neuroscience applications in order to obtain insights into the functionality of the brain. Finally, we consider the application of material informatics, learning from materials data. Here, we deploy regression techniques on materials data to predict physical properties of materials.

Matrix balancing is a preprocessing step in linear algebra computations such as the computation of eigenvalues of a matrix. Such computations are known to be numerically unstable if the matrix is unbalanced, that is the L2 norm of some rows and their corresponding columns are different by orders of magnitude. Given an unbalanced matrix A, the goal of matrix balancing is to find an invertible diagonal matrix D such that DAD^{−1} is balanced or approximately balanced in the sense that every row and its corresponding column have the same norm. In thesis, we study a classic iterative algorithm for matrix balancing due to Osborne (1960). The original algorithm was proposed for balancing rows and columns in the L2 norm, and it works by iterating over balancing a row-column pair in fixed round-robin order. Variants of the algorithm for other norms have been heavily studied and are implemented as standard preconditioners in many numerical linear algebra packages. Despite the popularity of Osborne’s algorithm in practice and extensive research on it there were no polynomial-time upper bound on the running time of this algorithm to explain the excellent performance of this algorithm in practice. Recently (in 2015), Schulman and Sinclair, in a first result of its kind for any norm, analyzed the rate of convergence of a variant of Osborne’s algorithm that uses the L∞ norm and a different order of choosing row-column pairs. In this thesis we study matrix balancing in the L1 norm and other Lp norms. We consider two notions of approximately balancing matrices and refer to them as ε-balancing and strict ε-balancing. As the names suggest strict ε-balancing implies ε-balancing. These notions will be defined in the body of the thesis.We prove upper bounds on the convergence rate of Osborne’s algorithm and some of its vari- ants. We prove fast convergence of different variants of the algorithm to an ε-balanced matrix, and propose a variant that converges to a strictly ε-balanced matrix in polynomial time. These results resolve a problem that has been open since Osborne proposed his algorithm in 1960. The following is a summary of our results for any real matrix A = (a_{ij})_{i,j}:1. We propose a simple greedy variant of Osborne’s algorithm and show that it converges to an ε-balanced matrix in K = O(min{ε^{−2} log w, ε^{−1}n^{3/2} log(w/ε)}) iterations that cost a total of O(m + Kn log n) arithmetic operations over O(n log(w/ε))-bit numbers. Here m is the number of non-zero entries of A, and w = \sum_{i,j} |a_{ij} |/a_{min} with a_{min}= min{|a_{ij} | : a_{ij} ̸= 0}.2. We show that the original round-robin implementation of Osborne’s algorithm converges to an ε-balanced matrix in O(ε^{−2}n^2 log w) iterations totaling O(ε^{−2}mn log w) arithmetic operations over O(n log(w/ε))-bit numbers.3. We devise a random implementation of the iteration and show that it converges to an ε-balanced matrix in O(ε^{−2} log w) iterations using O(m + ε^{−2}n log w) arithmetic operations over O(log(wn/ε))-bit numbers.4. We propose a variant of Osborne’s algorithm and prove that it converges to a strictly ε-balanced matrix in O (ε^{−2}n^9 log(wn/ε) log w/ log n) iterations that requireO (ε^{−2}n^10 log(wn/ε) log w/ log n) arithmetic operations over O(n log w/ε)-bit numbers.5. We demonstrate a lower bound of Ω(1/\sqrt{ε}) on the convergence rate of any implementation of the iteration. Thus, the dependence of our upper bounds on 1/ε is in the right ballpark.All our results are proved for balancing in L1 norm, but we observe, through a known trivial reduction, that our results for L1 balancing apply to any Lp norm for all finite p, at the cost of increasing the number of iterations by only a factor of p (or p2 in some cases).

Many scientific simulation codes represent physical systems by storing dataat a set of discrete points. For example, a weather modeling software applicationmay store data representing temperature, wind velocity, air pressure, humidity,etc. at a set of points in the atmosphere over a portion of the Earth's surface.Similarly, Computational Fluid Dynamics (CFD) software stores data representingfluid velocity, pressure, etc. at discrete points in the domain of the particularproblem being solved. In many simulations, this set of discrete points isa structured grid; i.e., a finite n-dimensional regular lattice.One of the most important steps in many scientific simulations is the solutionof a system of linear equations, where the unknowns of the system correspondto data elements at each of the discrete points used to model the system.If the simulation is based on a structured grid, this linear system willoften have a special structure itself, and this structure may lead tomore efficient techniques for solving the system than can be used forgeneral sparse linear systems.Scientific simulations most often use an iterative technique such asthe Conjugate Gradient Method or GMRES for solving linear systems..It is well known that these iterativetechniques converge much more quickly if a preconditioner is used. The ILU,or Incomplete LU Factorization, preconditioner is a good choice but it isnot parallelizable is a straightforward way.This thesis examines techniques for parallelizing the application of the ILUpreconditioner to linear systems arising from scientific simulations on structured grids.Various techniques are tested and timing results are recorded for different types andsizes of linear systems on structured grids.

High dimensional data usually have intrinsic low rank representations. These low rank representations not only reveal the hidden structure of the data but also reduce the computational cost of data analysis. Therefore, finding low dimensional approximations of the data is an essential task in many data mining applications.Classical low dimensional approximations rely on two universal tools: the eigenvalue decomposition and the singular value decomposition. These two different but related decompositions are of high importance in a large number of areas in science and engineering. As a result, research in numerical linear algebra has been conducted to derive efficient algorithms for solving eigenvalue and singular value problems. Because available solvers for these problems are so well developed, they are often used as black boxes in data analysis.This thesis explores numerical linear algebra techniques and extends well-known methods for low rank approximations to solve new problems in data analysis. Specifically, we carefully analyze the trace ratio optimization and argue that solving this problem can be done efficiently. We also propose efficient algorithms for low rank matrix approximations with missing entries. We also reformulate and analyze classical problems from a different perspective. This reveals the connection between the proposed methods and traditional methods in numerical linear algebra. The performance of the proposed algorithms is established both theoretically and through extensive experiments in dimension reduction and collaborative filtering.

While the theoretical analysis of linear dynamical systems with finite state-spaces is a mature topic, in situations where the underlying model has a large number of dimensions, modelers must turn to computational tools to better visualize and analyze the dynamic behavior of interest. In these situations, we are confronted with the Curse of Dimensionality: computational and storage complexity grows exponentially in the number of dimensions. This doctoral project focuses on two main classes of large-scale linear systems which arise in system biology. The Chemical Master Equation (CME) is a Fokker-Planck equation which describes the evolution of the probability mass function of a countable state space Markov process. Each state of the CME is labelled with an ordered S-tuple corresponding to one configuration of a well-mixed chemical system, where S is the number of distinct chemical species of interest. Even in cases where one only considers a projection of the CME to a finite subset of the states, one still must contend with the Curse of Dimensionality: the computational complexity grows exponentially in the number of chemical species. This dissertation describes a computational methodology for efficient solution of the CME which, in the best cases, will scale linearly in the number of chemical species. The second main class of high-dimensional problems requiring computational tools are coupled linear reaction-diffusion equations. For this class of models, we focus primarily on the computation of certain high-dimensional matrices which describe in a quantitative sense the input-to-state and state-to-output relationships. We describe algorithms for extracting useful information stored in these matrices and use this information to efficiently compute both reduced order models and open-loop control laws for steering the full system. A key feature of this approach is that the method is completely simulation or experiment free, in fact, in our numerical experiments, the computation of a reduced model or open-loop control law is an order of magnitude faster on a laptop than simulation of the full system on a 32 core node of a high-performance cluster. In both projects, the enabling computational technology is the recently proposed Tensor Train (TT) structured low-parametric representation of high-dimensional data. The TT-format effectively exploits low-rank structure of the "unfolding matrices" for compression and computational efficiency. Formally, the computational complexity of basic TT arithmetics scale linearly in the number of dimensions, potentially circumventing the curse of dimensionality. To demonstrate the effectiveness of this approach, we performed numerous numerical experiments whose results are reported here.

This thesis targets the design of parallelizable algorithms and communication-efficient parallel schedules for numerical linear algebra as well as computations with higher-order tensors. Communication is a growing bottleneck in the execution of most algorithms on parallel computers, which manifests itself as data movement both through the network connecting different processors and through the memory hierarchy of each processor as well as synchronization between processors. We provide a rigorous theoretical model of communication and derive lower bounds as well as algorithms in this model. Our analysis concerns two broad areas of linear algebra and of tensor contractions. We demonstrate the practical quality of the new theoretically-improved algorithms by presenting results which show that our implementations outperform standard libraries and traditional algorithms. We model the costs associated with local computation, interprocessor communication and synchronization, as well as memory to cache data transfers of a parallel schedule based on the most expensive execution path in the schedule. We introduce a new technique for deriving lower bounds on tradeoffs between these costs and apply them to algorithms in both dense and sparse linear algebra as well as graph algorithms. These lower bounds are attained by what we refer to as 2.5D algorithms, which we give for matrix multiplication, Gaussian elimination, QR factorization, the symmetric eigenvalue problem, and the Floyd-Warshall all-pairs shortest-paths algorithm. 2.5D algorithms achieve lower interprocessor bandwidth cost by exploiting auxiliary memory. Algorithms employing this technique are well known for matrix multiplication, and have been derived in the BSP model for LU and QR factorization, as well as the Floyd-Warshall algorithm. We introduce alternate versions of LU and QR algorithms which have measurable performance improvements over their BSP counterparts, and we give the first evaluations of their performance. We also explore network-topology-aware mapping on torus networks for matrix multiplication and LU, showing how 2.5D algorithms can efficiently exploit collective communication, as well as introducing an adaptation of Cannon's matrix multiplication algorithm that is better suited for torus networks with dimension larger than two. For the symmetric eigenvalue problem, we give the first 2.5D algorithms, additionally solving challenges with memory-bandwidth efficiency that arise for this problem. We also give a new memory-bandwidth efficient algorithm for Krylov subspace methods (repeated multiplication of a vector by a sparse-matrix), which is motivated by the application of our lower bound techniques to this problem. The latter half of the thesis contains algorithms for higher-order tensors, in particular tensor contractions. The motivating application for this work is the family of coupled-cluster methods, which solve the many-body Schrödinger equation to provide a chemically-accurate model of the electronic structure of molecules and chemical reactions where electron correlation plays a significant role. The numerical computation of these methods is dominated in cost by contraction of antisymmetric tensors. We introduce Cyclops Tensor Framework, which provides an automated mechanism for network-topology-aware decomposition and redistribution of tensor data. It leverages 2.5D matrix multiplication to perform tensor contractions communication-efficiently. The framework is capable of exploiting symmetry and antisymmetry in tensors and utilizes a distributed packed-symmetric storage format. Finally, we consider a theoretically novel technique for exploiting tensor symmetry to lower the number of multiplications necessary to perform a contraction via computing some redundant terms that allow preservation of symmetry and then cancelling them out with low-order cost. We analyze the numerical stability and communication efficiency of this technique and give adaptations to antisymmetric and Hermitian matrices. This technique has promising potential for accelerating coupled-cluster methods both in terms of computation and communication cost, and additionally provides a potential improvement for BLAS routines on complex matrices.

Sparse matrix formats encode very large numerical matrices with relatively few nonzeros. They are typically implemented using imperative languages, with emphasis on low-level optimization. Such implementations are far removed from the conceptual organization of the sparse format, which obscures the representation invariants. They are further specialized by hand for particular applications, machines and workloads. Consequently, they are difficult to write, maintain, and verify.In this dissertation we present LL, a small functional language suitable for implementing operations on sparse matrices. LL supports nested list and pair types, which naturally represent compressed matrices. It provides a few built-in operators and has a unique compositional dataflow model. As a result, LL programs are often direct analogs of high-level dataflow diagrams. This correspondence is useful when implementing sparse format construction and algebraic kernels such as sparse matrix-vector multiplication (SpMV). Despite its limited expressiveness, LL can specify interesting operations on real-world sparse formats.Next, we present a full functional verifier for sparse codes written in LL. We use a higher-order logic theorem prover, deploying standard simplification and introduction techniques. We use novel parametric predicates for tracking relationships between mathematical objects and their concrete representation. A simple heuristic tactic is used for automating proofs of different sparse formats. A qualitative analysis shows that our rule base exhibits good reusability and that the system easily extends to hierarchical formats.Finally, we describe a compiler for LL programs that generates efficient, parallel C code. We systematically map high-level nested datatypes onto compact and efficient low-level types. This facilitates a straightforward, syntax-directed translation of code. Local optimizations improve the performance of nested loops and eliminate redundant memory accesses. A coarse-grained heuristic loop parallelization scheme produces data-parallel code for execution on shared-memory multicore machines. Empirical evaluation shows that the sequential throughput of generated SpMV code is within 0.94--1.09 of a handwritten implementation. Parallel execution is 1.82--1.98 faster on two cores and 2.44--3.59 faster on four cores.

Data analysis is a process of inspecting and obtaining useful information from the data, with the goal of knowledge and scientific discovery. It brings together several disciplines in mathematics and computer science, including statistics, machine learning, database, data mining, and pattern recognition, to name just a few. A typical challenge with the current era of information technology is the availability of large volumes of data, together with ``the curse of dimensionality''. From the computational point of view, such a challenge urges efficient algorithms that can scale with the size and the dimension of the data. Numerical linear algebra lays a solid foundation for this task via its rich theory and elegant techniques. There are a large amount of examples which show that numerical linear algebra consists of a crucial ingredient in the process of data analysis. In this thesis, we elaborate on the above viewpoint via four problems, all of which have significant real-world applications. We propose efficient algorithms based on matrix techniques for solving each problem, with guaranteed low computational costs and high quality results. In the first scenario, a set of so called Lanczos vectors are used as an alternative to the principal eigenvectors/singular vectors in some processes of reducing the dimension of the data. The Lanczos vectors can be computed inexpensively, and they can be used to preserve the latent information of the data, resulting in a quality as good as by using eigenvectors/singular vectors. In the second scenario, we consider the construction of a nearest-neighbors graph. Under the framework of divide and conquer and via the use of the Lanczos procedure, two algorithms are designed with sub-quadratic (and close to linear) costs, way more efficient than existing practical algorithms when the data at hand are of very high dimension. In the third scenario, a matrix blocking algorithm for reordering and finding dense diagonal blocks of a sparse matrix is adapted, to identify dense subgraphs of a sparse graph, with broad applications in community detection for social, biological and information networks. Finally, in the fourth scenario, we visit the classical problem of sampling a very high dimensional Gaussian distribution in statistical data analysis. A technique of computing a function of a matrix times a vector is developed, to remedy traditional techniques (such as via the Cholesky factorization of the covariance matrix) that are limited to mega-dimensions in practice. The new technique has a potential for sampling Guassian distributions in tera/peta-dimensions, which is typically required for large-scale simulations and uncertainty quantifications./