A New Walk on Equations Monte Carlo Method for Linear Algebraic Problems

International audienceA new Walk on Equations (WE) Monte Carlo algorithm for Linear Algebra (LA) problem is proposed and studied. This algorithm relies on a non-discounted sum of an absorbed random walk. It can be applied for either real or complex matrices. Several techniques like simultaneous scoring or the sequential Monte Carlo method are applied to improve the basic algorithm. Numerical tests are performed on examples with matrices of different size and on systems coming from various applications. Comparisons with standard deterministic or Monte Carlo algorithms are also done

Asynchronous Approximation of a Single Component of the Solution to a Linear System

We present a distributed asynchronous algorithm for approximating a single component of the solution to a system of linear equations $Ax = b$, where $A$ is a positive definite real matrix, and $b \in \mathbb{R}^n$. This is equivalent to solving for $x_i$ in $x = Gx + z$ for some $G$ and $z$ such that the spectral radius of $G$ is less than 1. Our algorithm relies on the Neumann series characterization of the component $x_i$, and is based on residual updates. We analyze our algorithm within the context of a cloud computation model, in which the computation is split into small update tasks performed by small processors with shared access to a distributed file system. We prove a robust asymptotic convergence result when the spectral radius $\rho(|G|) < 1$, regardless of the precise order and frequency in which the update tasks are performed. We provide convergence rate bounds which depend on the order of update tasks performed, analyzing both deterministic update rules via counting weighted random walks, as well as probabilistic update rules via concentration bounds. The probabilistic analysis requires analyzing the product of random matrices which are drawn from distributions that are time and path dependent. We specifically consider the setting where $n$ is large, yet $G$ is sparse, e.g., each row has at most $d$ nonzero entries. This is motivated by applications in which $G$ is derived from the edge structure of an underlying graph. Our results prove that if the local neighborhood of the graph does not grow too quickly as a function of $n$, our algorithm can provide significant reduction in computation cost as opposed to any algorithm which computes the global solution vector $x$. Our algorithm obtains an $\epsilon \|x\|_2$ additive approximation for $x_i$ in constant time with respect to the size of the matrix when the maximum row sparsity $d = O(1)$ and $1/(1-\|G\|_2) = O(1)$

A Monte Carlo method for computing the action of a matrix exponential on a vector

A Monte Carlo method for computing the action of a matrix exponential for a certain class of matrices on a vector is proposed. The method is based on generating random paths, which evolve through the indices of the matrix, governed by a given continuous-time Markov chain. The vector solution is computed probabilistically by averaging over a suitable multiplicative functional. This representation extends the existing linear algebra Monte Carlo-based methods, and was used in practice to develop an efficient algorithm capable of computing both, a single entry or the full vector solution. Finally, several relevant benchmarks were executed to assess the performance of the algorithm. A comparison with the results obtained with a Krylov-based method shows the remarkable performance of the algorithm for solving large-scale problems.info:eu-repo/semantics/acceptedVersio

A distributed Monte Carlo based linear algebra solver applied to the analysis of large complex networks

Methods based on Monte Carlo for solving linear systems have some interesting properties which make them, in many instances, preferable to classic methods. Namely, these statistical methods allow the computation of individual entries of the output, hence being able to handle problems where the size of the resulting matrix would be too large. In this paper, we propose a distributed linear algebra solver based on Monte Carlo. The proposed method is based on an algorithm that uses random walks over the system’s matrix to calculate powers of this matrix, which can then be used to compute a given matrix function. Distributing the matrix over several nodes enables the handling of even larger problem instances, however it entails a communication penalty as walks may need to jump between computational nodes. We have studied different buffering strategies and provide a solution that minimizes this overhead and maximizes performance. We used our method to compute metrics of complex networks, such as node centrality and resolvent Estrada index. We present results that demonstrate the excellent scalability of our distributed implementation on very large networks, effectively providing a solution to previously unreachable problem instances.info:eu-repo/semantics/acceptedVersio

A probabilistic linear solver based on a multilevel Monte Carlo Method

We describe a new Monte Carlo method based on a multilevel method for computing the action of the resolvent matrix over a vector. The method is based on the numerical evaluation of the Laplace transform of the matrix exponential, which is computed efficiently using a multilevel Monte Carlo method. Essentially, it requires generating suitable random paths which evolve through the indices of the matrix according to the probability law of a continuous-time Markov chain governed by the associated Laplacian matrix. The convergence of the proposed multilevel method has been discussed, and several numerical examples were run to test the performance of the algorithm. These examples concern the computation of some metrics of interest in the analysis of complex networks, and the numerical solution of a boundary-value problem for an elliptic partial differential equation. In addition, the algorithm was conveniently parallelized, and the scalability analyzed and compared with the results of other existing Monte Carlo method for solving linear algebra systems.info:eu-repo/semantics/acceptedVersio

A highly parallel algorithm for computing the action of a matrix exponential on a vector based on a multilevel Monte Carlo method

A novel algorithm for computing the action of a matrix exponential over a vector is proposed. The algorithm is based on a multilevel Monte Carlo method, and the vector solution is computed probabilistically generating suitable random paths which evolve through the indices of the matrix according to a suitable probability law. The computational complexity is proved in this paper to be significantly better than the classical Monte Carlo method, which allows the computation of much more accurate solutions. Furthermore, the positive features of the algorithm in terms of parallelism were exploited in practice to develop a highly scalable implementation capable of solving some test problems very efficiently using high performance supercomputers equipped with a large number of cores. For the specific case of shared memory architectures the performance of the algorithm was compared with the results obtained using an available Krylov-based algorithm, outperforming the latter in all benchmarks analyzed so far.info:eu-repo/semantics/acceptedVersio

A Fast Monte Carlo algorithm for evaluating matrix functions with application in complex networks

We propose a novel stochastic algorithm that randomly samples entire rows and columns of the matrix as a way to approximate an arbitrary matrix function. This contrasts with the "classical" Monte Carlo method which only works with one entry at a time, resulting in a significant better convergence rate than the "classical" approach. To assess the applicability of our method, we compute the subgraph centrality and total communicability of several large networks. In all benchmarks analyzed so far, the performance of our method was significantly superior to the competition, being able to scale up to 64 CPU cores with a remarkable efficiency.Comment: Submitted to the Journal of Scientific Computin

State estimation and trajectory tracking control for a nonlinear and multivariable bioethanol production system

In this paper a controller is proposed based on linear algebra for a fed-batch bioethanol production process. It involves fnding feed rate profles (control actions obtained as a solution of a linear equations system) in order to make the system follow predefned concentration profles. A neural network states estimation is designed in order to know those variables that cannot be measured. The controller is tuned using a Monte Carlo experiment for which a cost function that penalizes tracking errors is defned. Moreover, several tests (adding parametric uncertainty and perturbations in the control action) are carried out so as to evaluate the controller performance. A comparison with another controller is made. The demonstration of the error convergence, as well as the stability analysis of the neural network, are included.Fil: Fernández, Maria Cecilia. Universidad Nacional de San Juan. Facultad de Ingeniería. Instituto de Ingeniería Química; ArgentinaFil: Pantano, Maria Nadia. Universidad Nacional de San Juan. Facultad de Ingeniería. Instituto de Ingeniería Química; Argentina. Consejo Nacional de Investigaciones Científicas y Técnicas; ArgentinaFil: Rossomando, Francisco Guido. Consejo Nacional de Investigaciones Científicas y Técnicas. Centro Científico Tecnológico Conicet - San Juan. Instituto de Automática. Universidad Nacional de San Juan. Facultad de Ingeniería. Instituto de Automática; ArgentinaFil: Ortiz, Oscar Alberto. Universidad Nacional de San Juan. Facultad de Ingeniería. Instituto de Ingeniería Química; ArgentinaFil: Scaglia, Gustavo Juan Eduardo. Universidad Nacional de San Juan. Facultad de Ingeniería. Instituto de Ingeniería Química; Argentina. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentin

A Stochastic Method for Solving Time-Fractional Differential Equations

We present a stochastic method for efficiently computing the solution of time-fractional partial differential equations (fPDEs) that model anomalous diffusion problems of the subdiffusive type. After discretizing the fPDE in space, the ensuing system of fractional linear equations is solved resorting to a Monte Carlo evaluation of the corresponding Mittag-Leffler matrix function. This is accomplished through the approximation of the expected value of a suitable multiplicative functional of a stochastic process, which consists of a Markov chain whose sojourn times in every state are Mittag-Leffler distributed. The resulting algorithm is able to calculate the solution at conveniently chosen points in the domain with high efficiency. In addition, we present how to generalize this algorithm in order to compute the complete solution. For several large-scale numerical problems, our method showed remarkable performance in both shared-memory and distributed-memory systems, achieving nearly perfect scalability up to 16,384 CPU cores.Comment: Submitted to the Journal of Computational Physic

A Study of Biased and Unbiased Stochastic Algorithms for Solving Integral Equations

In this paper we propose and analyse a new unbiased stochastic method for solving a class of integral equations, namely the second kind Fredholm integral equations. We study and compare three possible approaches to compute linear functionals of the integral under consideration: i) biased Monte Carlo method based on evaluation of truncated Liouville-Neumann series, ii) transformation of this problem into the problem of computing a finite number of integrals, and iii) unbiased stochastic approach. Five Monte Carlo algorithms for numerical integration have been applied for approach (ii). Error balancing of both stochastic and systematic errors has been discussed and applied during the numerical implementation of the biased algorithms. Extensive numerical experiments have been performed to support the theoretical studies regarding the convergence rate of Monte Carlo methods for numerical integration done in our previous studies. We compare the results obtained by some of the best biased stochastic approaches with the results obtained by the proposed unbiased approach. Conclusions about the applicability and efficiency of the algorithms have been drawn