A Note on PA = LU in Javascript

• Download lu-decomposition.zip - 14.5 KB

Introduction

This is a note on LU decomposition in Javascript.

Background

Through out the human history, linear simultaneous equations had and will continue to have such a significant influence to our day to day life. It is the core component of the solutions to a large section of more complicated numerical computation problems. The solution to the linear simultaneous equations is probably the best studied subject in the human history.

• The famous Carl Friedrich Gauss - Gauss studied this subject and the "Gaussian elimination" is named after him;
• The famous Isaac Newton - Newton also studied this subject and the Cambridge University published his notes as Arithmetica Universalis in 1707. Because Newton is more than 100 years older than Gauss, he obviously touched this subject earlier;
• The not very famous Liu Hui - Liu also studied this subject and wrote a book "The Nine Chapters on the Mathematical Art". Because Liu is more than 1000 years older than Newton, he obviously touched this subject earlier. Since part of the chapters were written in approximately 150 BC, humans may have studied the solution to linear simultaneous equations before Jesus Christ was born;
• The famous Alan Turing - Turing studied this subject with the context of the digital computers. He invented the "LU decomposition" based on the well known "Gaussian elimination", which is well suited to solve the linear simultaneous equations by computers.

Solving linear simultaneous equations is such a important subject and is also a relatively easy one. It is a good chance to have some daily exercise to prevent the Alzheimer's Disease. So let us implement the famous algorithm PA = LU in Javascript and have some fun with it.

The attached is a Java Maven project. But you do not need to have Java to run it. The example is intended to run in a web browser, so the important files are the Javascript and the HTML files. Since the "index.html" file referenced a CDN, you will need a computer that has internet access to load it.

The Algorithm

The standard implementation of the "Gaussian elimination" in the context of the digital computers is the "LU decomposition" introduced by Alan Turing. Typically we need the pivoting operations in the LU decomposition. But for simplicity reasons, let first take a look at the steps of a basic LU decomposition.

The Basic LU Decomposition

For a given matrix A, the goal of the LU decomposition is to find a lower diagonal matrix L and an upper diagonal matrix U, such that A = LU.

In order to find the L and U, we start by setting L as an idenity matrix and U = A.

As the first step to make the U matrix upper diagonal, we need to eliminate all the first column entries under the first row.

It is not difficult to verify that A = LU holds after the operations. The second step of the LU decomposition is to eliminate the entries in the second column under the second row.

We can continue the process to eliminate all the entries under the diagonal entries to complete the LU decomposition for A = LU.

PA = LU

In many cases, we are not so lucky that all the diagonal entries are non-zero at all the steps. Even when the matrix A is non-singular, zero diagonal entries can exist. In such cases, we need to swap the rows in the U matrix to create a non-zero diagonal entry. The process of swapping the rows is called "Pivoting" or "permutation".

In order to maintain the equality after swapping the rows, a permutation matrix P is introduced. At the 0_th step of the LU decomposition, the P matrix is initialized as an identity matrix.

• Permutation is needed when the diagonal entry is 0, because it is the divisor in the process of the LU decomposition;
• When the diagonal entry is non-zero, permutation is needed for numerical stability reasons. This is because the numerical error can be significantly larger if the divisor is too small.

In the example above, we will scan all the entries from u_33⁽²⁾ to u_n3⁽²⁾ to find the entry with the largest absolute value, which is u_p3⁽²⁾ in this case. We will swap the row 3 and row p in the U matrix to move u_p3⁽²⁾ to the diagonal. In order to maintain PA = LU, we need to mutate the P and L matrices at the same time. In order to see how we need to mutate the P and L matrices, let us sub-divide the L and U matrices into block matrices.

From the above block matrix operations, we can see that swapping the row 3 and row p in the U matrix only affect the lower part of the matrix product LU. With this in mind, we can use the following fact to figure out how we can mutate the P and L matrices.

If you take a little detailed look into the block matrices, it is not difficult to find that the operations to maintain PA = LU and to move the u_p3⁽²⁾ to the diagonal are the following.

• Swap the row 3 and row p in the U matrix for columns 3 to n;
• Swap the row 3 and row p in the L matrix for columns 1 to 2;
• Swap the row 3 and row p in the P matrix.

Regardless if the diagonal entry is zero, pivoting is typically needed for better numerical stability for every elimination step of the LU decomposition.

The Pseudo Code PA = LU

In summary, the algorithm for LU decomposition with partial pivoting PA = LU can be described by the following pseudo code.

The Javascript Implementation

The Matrix Representation in Javascript

In my implementation of the algorithms, a matrix in Javascript is implemented by an array of arrays.

• Each entry in the first layer of the array represents a row in the matrix;
• Each entry in a second layer of the array represents the entry of the corresponding column in the row.

Because I am using Javascript, all the indices of the array elements are 0 based.

The Utilities

In order to make the work a little easier, I implemented an utility object that I will be using extensively in the algorithm implementation. The utility object exposes a few methods. The purpose of each method is easily identified by it name.

// Utility
let MatrixUtil = function() {
    let Identity =  function(dim) {
        let I = [];
    
        let i, j;
        for (i = 0; i < dim; i++) {
            let r = [];
            
            for (j = 0; j < dim; j++) {
                r.push((i == j)? 1: 0);
            }
            
            I.push(r)
        }
        
        return I;
    };
    
    let Clone = function(M) {
        let C = [];
        let i, j;
        
        let rc = M.length;
        for (i = 0; i < rc; i++) {
            let m = M[i];
            let r = [];
            
            let cc = m.length;
            for (j = 0; j < cc; j++) {
                r.push(m[j]);
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    let Transpose = function(M) {
        let C = [];
        let i, j;
        
        let rc = M.length;
        let cc = M[0].length;
        
        for (i = 0; i < cc; i++) {
            let r = [];
            
            for (j = 0; j < rc; j++) {
                r.push(M[j][i]);
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    let Sum = function(A, B) {
        let C = [];
        let i, j;
        
        let rc = A.length;
        let cc = A[0].length;
    
        for (i = 0; i < rc; i++) {
            let r = [];
            
            for (j = 0; j < cc; j++) {
                r.push(A[i][j] + B[i][j]);
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    let Subtract =  function(A, B) {
        let C = [];
        let i, j;
        
        let rc = A.length;
        let cc = A[0].length;
        for (i = 0; i < rc; i++) {
            let r = [];
            
            for (j = 0; j < cc; j++) {
                r.push(A[i][j] - B[i][j]);
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    let Multiply =  function(A, B) {
        let C = [];
        
        let rCount = A.length;
        let cCount = B[0].length;
        
        let i, j, k;
        for (i = 0; i < rCount; i++) {
            let r = [];
            let ra = A[i];
            
            for (j = 0; j < cCount; j++) {
                
                let cell = 0;
                for (k = 0; k < rCount; k++) {
                    cell = cell + ra[k] * B[k][j];
                }
                
                r.push(cell);
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    let Round = function(M, decimal) {
        let C = [];
        
        let rc = M.length;
        let cc = M[0].length;
        
        let i, j;
        for (i = 0; i < rc; i++) {
            let r = [];
            let rm = M[i];
            
            for (j = 0; j < cc; j++) {
                r.push(Number(rm[j].toFixed(decimal)));
            }
            
            C.push(r);
        }
        
        return C;
    };
    
    return {
        Identity: Identity,
        Clone: Clone,
        Transpose: Transpose,
        Sum: Sum,
        Subtract: Subtract,
        Multiply: Multiply,
        Round: Round
    };
}();

The PA = LU

The LU decomposition algorithm is implemented in the constructor function "LU".

let LU = function() {
    let self = this;
    
    let getMaxEntry = function(cIdex) {
        let U = self.U;
        let maxIdex = cIdex, max = Math.abs(U[cIdex][cIdex]);
        let dim = A.length, i;
        
        for (i = cIdex + 1; i < dim; i++) {
            let next = Math.abs(U[i][cIdex]);
            if (next > max) {
                max = next;
                maxIdex = i;
            }
        }
        
        return maxIdex;
    };
    
    let pivot = function(p, n) {
        let dim = self.A.length;
        let U = self.U, L = self.L, P = self.P;
        
        let temp, i;
        // U
        for (i = p; i < dim; i++) {
            temp = U[p][i]; U[p][i] = U[n][i]; U[n][i] = temp;
        }
        
        // L
        for (i = 0; i < p; i++) {
            temp = L[p][i]; L[p][i] = L[n][i]; L[n][i] = temp;
        }
        
        // P
        temp = P[p]; P[p] = P[n]; P[n] = temp;
    };
    
    let eliminate = function(p) {
        let dim = self.A.length;
        let U = self.U, L = self.L;
        let i, j;
        
        for (i = p + 1; i < dim; i++) {
            let l = U[i][p] / U[p][p];
            L[i][p] = l;
            
            for (j = p; j < dim; j++) {
                U[i][j] = U[i][j] - l * U[p][j];
            }
        }
        
    };
    
    let setA = function(A) {
        // A is a square matrix
        let dim = A.length;
        
        self.A = A;
        self.U = MatrixUtil.Clone(A);
        self.L = MatrixUtil.Identity(dim);
        self.P = MatrixUtil.Identity(dim);
        
        return self;
    };
    
    let PLU = function() {
        let A = self.A;
        let dim = A.length;
        
        let i, j, k;
        for (i = 0; i < dim - 1; i++) {
            // Find the max entry
            let maxIdex = getMaxEntry(i);
            
            // Pivoting
            if (i != maxIdex) {
                pivot(i, maxIdex);
            }
            
            // Eliminate
            eliminate(i);
        }
    };
    
    let solve = function(B) {
        let dim = self.A.length;
        let P = self.P;
        let L = self.L;
        let U = self.U;
        
        let PB = MatrixUtil.Multiply(P, B);
        
        let i, j;
        
        //LUX = PB
        //UX = Y ==> LY = PB
        let Y = [];
        for (i = 0; i < dim; i++) {
            let l = L[i];
            
            let temp = 0;
            for (j = 0; j < i; j++) {
                temp = temp + l[j] * Y[j][0];
            }
            
            Y.push([PB[i][0] - temp]);
        }
        
        // UX = Y
        let X = [];
        for (i = dim - 1; i >= 0; i--) {
            let u = U[i];
            
            let start = i + 1;
            let temp = 0;
            for (j = start; j < dim; j++) {
                temp = temp + u[j] * X[j - start];
            }
            
            X.unshift([(Y[i][0] - temp)/u[i]]);
        }
        
        return X;
    };
    
    // Public methods
    self.SetA = setA;
    self.PLU = PLU;
    self.Solve = solve;
};

This constructor function exposes three method.

• The "SetA()" - it allows you to set the A matrix. It also initiates the P, L, and U matrices;
• The "PLU()" - it performs the LU decomposition on the A matrix;
• The "Solve()" - It takes a column vector as B and solves the equation AX = B. You need to make sure that it is called after the "PLU()" method is called.

Because we are using Javascript, you have access to the P, L, and U matrices through the instance of the LU object after the "PLU()" method is called for the A matrix set by the "SetA()" method.

The Example

The Example "index.html"

In order to give a try to my implementation of the LU decomposition, I created a small example program.

<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>lu-decomposition-example</title>
<link rel="stylesheet" type="text/css" href="styles/app.css">
    
<script type="text/javascript" src="script/MatrixUtil.js"></script>
<script type="text/javascript" src="script/lu-decomposition.js"></script>
    
<script src="https://cdnjs.cloudflare.com/ajax/libs/react/15.6.1/react.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/react/15.6.1/react-dom.min.js"></script>
    
<script type="text/javascript">
    // The example equition AX = B
    let A = function() {
        let A = [];
        A.push([2,1,1,0]);
        A.push([4,3,3,1]);
        A.push([8,7,9,5]);
        A.push([6,7,9,8]);
        
        return A;
    }();
    
    let B = function() {
        return [[5], [8], [1], [7]];        
    }();
    
    // Matrix DOM
    let React_Matrix = React.createClass({
        render: function() {
            let M = MatrixUtil.Round(this.props.M, 2);
            
            return React.DOM.table({className: 'matrix bb'},
                M.map(function(row, i) {
                    return (React.DOM.tr({key: i},
                        row.map(function(e, j) {
                            let p = {key: j};
                            let s = {style: {'padding-right': '1px',
                                visibility: (e < 0)? 'visible' : 'hidden'}};
                            
                            return React.DOM.td(p,
                                    React.DOM.span(s, '-'),
                                    React.DOM.span(null, Math.abs(e)));
                        })
                    ));
                })
            );
        }
    });
    
    let React_Main = React.createClass({
        getInitialState: function() {
            let lu = new LU().SetA(A);
    
            return { lu: lu, A: lu.A, B: B, P: lu.P, L: lu.L, U: lu.U };
        },
        LUandSolve: function() {
            let state = this.state;
            let lu = this.state.lu;
            lu.PLU();
            
            let X = lu.Solve(state.B);
            
            let R = MatrixUtil.Subtract(B, MatrixUtil.Multiply(state.A, X));
            let RT = MatrixUtil.Transpose(R);
            let Residual = MatrixUtil.Multiply(RT, R)[0][0];
            this.setState({ P: lu.P, L: lu.L, U: lu.U, X: X, Residual:  Residual});
        },
        Reset: function() {
            let lu = this.state.lu;
            lu.SetA(A);
    
            this.setState({ P: lu.P, L: lu.L, U: lu.U, X: null, Residual: null });
        },
        render: function() {
            let state = this.state;
            
            let e = [];
            e.push(React.DOM.div(null,
                    React.DOM.button({ onClick: this.LUandSolve, className: 'btn' },
                            'LU & Solve Equition'),
                    React.DOM.button({ onClick: this.Reset, className: 'btn' }, 'Reset')));
            
            let X = state.X;
            let Residual = state.Residual;
            
            e.push(React.DOM.div({className: 'y-center'},
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'A'),
                            React.createElement(React_Matrix, {M: state.A})),
                            
                    (X != null)? React.DOM.div(null, React.DOM.div({className: 'semi-bold'},
                            'X'), React.createElement(React_Matrix, {M: state.X}))
                            : React.DOM.span({className: 'semi-bold',
                                style: {padding: '0px 5px'}}, 'X'),
                    
                    React.DOM.span({style: {padding: '0px 5px'}}, '='),
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'B'),
                            React.createElement(React_Matrix, {M: state.B}))));
            
            e.push(React.DOM.div(null, (Residual == null)? null
                    :React.DOM.span({style: {padding: '0px 5px'}},
                            'Residual = ' + Residual)
            ));
            
            e.push(React.DOM.div({className: 'y-center'},
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'P'),
                        React.createElement(React_Matrix, {M: state.P})),
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'A'),
                        React.createElement(React_Matrix, {M: state.A})),
                    React.DOM.span({style: {padding: '0px 5px'}}, '='),
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'L'),
                        React.createElement(React_Matrix, {M: state.L})),
                    React.DOM.div(null, React.DOM.div({className: 'semi-bold'}, 'U'),
                        React.createElement(React_Matrix, {M: state.U}))));
    
            
            return React.DOM.div(null, e);
        }
    });
    
    
        
    window.onload = function() {
        ReactDOM.render(React.createElement(React_Main),
                document.getElementById('example'));
    };
    
</script>
</head>
<body>
    
<div id="example"></div>
    
</body>
</html>

This example uses React to display the matrices. If you are not familiar with React, you can take look at my early notes or go to the official React web pages.

• The A and B matrices are initiated at the beginning of the program;
• Two buttons are added to the web page. The "LU & Solve Equition" button will perform the LU decomposition and solve the equation AX = B. The "Reset" button will reset the result and let you to try it again;
• All the results, including P, L, U, and the solution of the equation is displayed through the React on the browser.

Run the Example

If you are not using a very old browser, you should be able to load the example and try it. Because the web page references a CDN for the React libraries, your computer needs to have internet access.

When the example is loaded, All the matrices are initiated. If you click the "LU & Solve Equition", an LU decomposition is performed on the A matrix and the AX = B equation is solved.

You can see the solution for X and the residual for the solution of the equation. You can also click the "Reset" button to clear the result and try it again.

Points of Interest

• This is a note on LU decomposition in Javascript;
• The method to solve linear simultaneous equations is such a important subject. The LU decomposition is such a great tool that the mathematicians give to us. Obviously this note is not the end of the story. We still have a lot of questions, such as how to deal with singular matrices and how to take advantage of the sparsity of the matrices. But these questions are out of the scope of this note now;
• I hope you like my postings and I hope this note can help you one way or the other.

History

First Revision - 8/26/2017.