Alexandrov, Vassil Nikolov

Cadenas Medina, Jose Oswaldo

Megson, Graham M

Plaks, T P

English

Central Archive at the University of Reading

A linear algebra processor using Monte Carlo methodsT. P. Plaks*, G. M. Megson**, O. Cadenas** and V. N. Alexandrov*** London South Bank University** The University of Reading, UKIntroductionMany problems in control engineering, scientiﬁc computations, military and space applications can bereduced to the solution of basic Linear Algebra problems, for solving systems of linear equations: AX = B.For example many problems arising from partial diﬀerential equations where local eﬀects or ordering ofﬁnite diﬀerence/ﬁnite elements produce recurrences with limited spatial and temporal extension. Becausethe problems are often large and involve many zero elements there is great need to speed up computationsor exploit structure of the matrices [13, 8]. In the former case to allow larger problems to be solved ormore iterations of a problem with varying parameters to be solved in the same time. In the sparse caseboth computation time and hardware involved can be optimized by removing redundancies from standarddense matrix approaches to the problem.There are many ways to solve the system, AX = B [1, 7, 13], among which is to ﬁnd ﬁrst the inversionA−1 and then calculate X = A−1B. The inversion of matrices with Banded structure [1, 7] is the classicproblem in Linear Algebra. In this paper we consider of hardware implementation of these algorithmsand, therefore, we will use the systolic or regular array approach for implementation of the proposedmethods. The best known arrays for matrix inversion are based on the so-called Algebraic Path problem[6, 11, 12] and formulation of the Gauss Jordan method [3]. These arrays generally require 4n or 5n stepsusing n2 processors and do not exploit any pivoting. Unfortunately the arrays cannot be easily adaptedto the banded or sparse cases and require the same run-time and processors as for dense systems evenfor triangular matrices.In this paper we introduce a Linear Algebra Processor based on Monte Carlo method for matrixinversion [4, 5, 14]. The processor is based on a regular array approach and is implemented on an FPGAplatform. The most important part of the processor is the Matrix inversion regular processor array basedon Monte Carlo method. The Monte Carlo array uses a modiﬁcation of the technique proposed in [9]to exploit the zero structure in band systems [10]. This analysis results both in decreases in run-timeand array latency but also in drastically reduced hardware compared with previous Monte-Carlo arrayfor dense matrices and compared to previous solvers. In particular we consider arrays with O(TNw)to O(nTNw) elements where T < √n is the length of a Markov chain, N is the number of chains, nis the order of the matrix and w is the band width. The run-time ranges from O(n2 + T + N + w) toO(n + T + N + w) when the design is laid out in higher dimensions. Consequently we produce arrayswith computation time t < 3n+ w steps where a step is the cost of a inner product.A part of the Monte-Carlo based matrix inversion regular array has been implemented in Virtex-II to invert a matrix of size 64 × 64. The FPGA circuit implementation runs at a frequency of over200 MHz. It is observed that a well structured and regular architecture can balance the computationalresources of modern FPGA while maintaining high performance. The architecture is not limited byFPGA computational resources to invert matrices of moderate size, but the practical implementation ofthe architecture for matrices of large sizes is limited by the number of pin-out count of the design [2].Monte Carlo methodDense matricesTo ﬁnd the inverse C = A−1 of a dense matrix B of order n we compute D = I − A (I is the identitymatrix), and then use the Monte Carlo method described below [9]:crr′ ≈ 1NN∑s=1 ∑(j|kj=r′)Wjs, Wj = Wj−1dkj−1kjpkj−1kj, (1)where crr′ is the element of C. The sequence of indices kj , 0≤j≤T , form the states of a Markovchain with transition probabilities pkj−1kj from state kj−1 to kj . Transitions will be denoted as tuples< kj−1, kj >. Thus, a transition determines the element of D, dkj−1,kj . For each j-th state the Wj isdeﬁned recurrently. A Markov chain starts with a row index, i.e. for computing an element crr′, the chainstarts at k0 = r, and only these Wj are included into sum of crr′ for which kj = r′. To compute one rowof inverse, N Markov chains (trials) are used over which the averages are computed to get the result. AMarkov chain is denoted by index s, 1≤s≤N . Notice, that N and T are the parameters of a stochasticprocess: the mathematical expectations of the number and length of Markov chains respectively. Thus,N and T do not depend on the matrix size, but depend only on the required precision of computationand the norm of matrix [9]. In the following, we introduce matrix E such that ep,q = dp,q/pp,q.Banded and sparse matricesThe Monte Carlo method, because of the probabilistic nature, produces the possibility to cut oﬀ thezero parts of computations. Observe that the Monte Carlo methods operates with rows and columnsof the matrix not as objects with ﬁxed ordered structures but as sets of elements which are selectedfor computations in an arbitrary order. Thus, the number of elements in the sets does not aﬀect thecomputation structure and zero elements can be easily excluded from the set [10]. Suppose, the size ofmatrix is 1024× 1024, and there are 1% non-zero elements. Thus, there are about 10000 elements thatcan be rearranged as a new, dense matrix of size 100× 100.ConclusionIn this paper we present a linear algebra processor that is based on a novel, Monte Carlo method approachfor matrix inversion. This approach allows us to achieve signiﬁcant improvement in performance and toreduce the amount of required hardware. Also, the same approach is easily used for banded and sparsematrices reducing signiﬁcantly the time and hardware resources.References[1] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed Computation. Numerical Methods.Prentice/Hall, Int., 1989.[2] O. Cadenas, G. Megson, and T. Plaks. FPGA circuits for a Monte-Carlo based matrix inversionarchitecture. In T. P. Plaks and P. M. Athanas, editors, ERSA’03. The International Conference onEngineering of Reconfigurable Systems and Algorithms, June 23–26, 2003, Las Vegas, USA, page .CSREA Press, 2003.[3] M. Cosnard. Designing parallel algorithms for linearly connected processors and systolic arrays.Advances in Parallel Computing, 1:273–317, 1990.2[4] J. H. Curtiss. Monte Carlo methods for the iteration of linear operators. J. Math. Phys,32(4):209–232, 1954.[5] J. H. Curtiss. A theoretical comparision of the eﬃciencies of two classical methods and Monte Carlomethod for computing one component of the solution of a set of linear algebraic equations. In Proc.Symp. on Monte Carlo Methods, pages 191–233. John Wiley and Sons, 1956.[6] J. M. Delosme. A parallel algorithm for the algebraic path problem. In M. C. et al., editor, Paralleland Distributed Algorithms, pages 67–78. Bonas, France, 1988.[7] T. L. Freeman. Parallel Numerical Algorithms. Series in Computer Science. Prentice/Hall, Int.,1992.[8] G. H. Golub and P. van Dooren, editors. Numerical Linear Algebra, Digital Signal Processing andParallel Algorithms. NATO ASI Series F: Computer and Systems Sciences Vol. 70. Springer-Verlag,1991.[9] G. M. Megson, V. N. Aleksandrov, and I. T. Dimov. Systolic matrix inversion using a Monte Carlomethod. Parallel Algorithms and Applications, 3: 311–330, 1994.[10] T. P. Plaks, G. M. Megson, and V. N. Alexandrov. A family of eﬃcient regular arrays for bandand sparse matrix inversion using Monte Carlo methods. In M. Deville and R. Owens, editors,Proc. of the 16th IMACS World Congress 2000 on Scientific Computation, Applied Mathematicsand Simulation, Lausanne, Switzerland, August 21-25, 2000. (CD-ROM), 6 pages. IMACS, USA,2000.[11] S. Rajopadhye. An improved systolic algorithm for the algebraic path problem. INTEGRATION,pages 279–296, 1993.[12] T. Risset and Y. Robert. Synthesis of processor arrays for the algebraic path problem unifying oldresults and deriving new architectures. Parallel Processing Letters, 1(1):19–28, 1991.[13] H. A. van der Vorst and P. van Dooren, editors. Parallel Algorithms for Numerical Linear Algebra.Advances in Parallel Computing, 1. North-Holland, Amsterdam, 1990.[14] J. R. Westlake. A Handbook of Numerical Matrix Inversion and Solution of Linear Equations. JohnWiley and Sons, New York, 1968.3

A linear algebra processor using Monte Carlo methods

Abstract

Similar works

Full text

Available Versions

Central Archive at the University of Reading