Functional programming as an example of working with matrices from the theory of linear algebra

Introduction

The basis of working with matrices (in this article we will consider only two-dimensional matrices) is a powerful mathematical theory from the field of linear algebra. One definition or action follows from another; one function calls another. Therefore, for functional implementation of the functional of mathematical operations on matrices, functional languages ​​fit very well. In this article, we will look at specific examples in F # and give detailed comments on how this works. Since F # is part of the .NET family, the resulting functionality can be used without any problems in other imperative languages, for example C #.

Matrix definition and implementation in F #

Matrices are the basic and most important part of linear algebra. Matrices are often used in programming, for example, in 3D modeling or game development. Of course, the bicycle has long been invented and the necessary frameworks for working with matrices are already ready, and they can and should be used. This article does not aim at inventing a new framework, but shows the implementation of basic mathematical operations for working with matrices in a functional style using the F # programming language. As we review the material, we will turn to the mathematical theory of matrices and see how it can be implemented in code.

First, let's recall what a matrix is? Theory tells us the following

A rectangular table of numbers containing m rows and n columns is called a matrix of size mxn

Matrices are usually denoted by capital letters of the Latin alphabet and are written as

$$

$$A=\begin{pmatrix}a11&a12&...&a1n\\a21&a22&...&a2n\\...&...&...&...\\am1&am2&...&amn\end{pmatrix}$$

Or shortly

$$

$$A=(aij)$$

To work with matrices in F # we will create a record based on a two-dimensional table, which we will further expand with useful methods in order to perform the necessary mathematical operations on it.

type Matrix = { values: int[,] }
    with
        //    
    end

Add a helper method to initialize the recording with a two-dimensional array

static member ofArray2D (values: int [,]) = 
    { values = values }

The input argument to the function will be a two-dimensional array, and at its output, a record of type Matrix . Below is an example of record initialization.

let a = array2D [[1;0;2]
                 [3;1;0]]
let A = Matrix.ofArray2D a

Two matrices A = (aij) and B = (bij) are called equal if they coincide elementwise, i.e. aij = bij for all i = 1,2 ..., m and j = 1,2 ... n

To implement this rule, we will use the overridden operator == and add a couple of useful functions, which we will also need in the future.

static member sizes matrix =
    let rows = matrix.values.[*,0].Length
    let cols = matrix.values.[0,*].Length
    (rows, cols)

static member isEquallySized matrix1 matrix2 =
    let dim1 = Matrix.sizes matrix1
    let dim2 = Matrix.sizes matrix2
    (dim1 = dim2)

static member (==) (matrix1, matrix2) =
    if not (Matrix.isEquallySized matrix1 matrix2) then false
    else
        not (matrix1.values
               |> Array2D.mapi (fun x y v -> if matrix2.values.[x, y] <> v then false else true)
               |> Seq.cast<bool>
               |> Seq.contains false)

Let's take a closer look at the code above. As you can see, there are three functions. The first sizes function returns the dimension of the matrix as a tuple. Since we work only with rectangular matrices, to get the number of rows, we take a full slice of the first column and return its length.

let rows = matrix.values.[*,0].Length

Determining the number of columns works in a similar way - we get a full slice of the first row and return its length.

The following isEquallySized function compares the dimension of two matrices and returns true if they are equal. To do this, it uses the ready-made sizes function and simply compares the results.

The == operator for elementwise comparing two matrices seems more complicated, but now you will see that it is also simple.

Before comparing two matrices, we compare their dimension. If they are not equal, then there is no further point in checking, since it is already clear that the matrices will not be equal.

if not (Matrix.isEquallySized matrix1 matrix2) then false

Further, on the basis of the original matrices matrix1 and matrix2, we form a new matrix filled with true or false , depending on whether the corresponding cells of both matrices coincide.

matrix1.values
|> Array2D.mapi (fun x y v -> if matrix2.values.[x, y] <> v then false else true

The Array2D.mapi function iterates over all elements of matrix1 and passes three parameters
x to the handler - row index
y - column index
v - cell
contents We compare the contents of cell v with the corresponding cell matrix2 and if they are equal, then write true , otherwise false .

If there is at least one cell with the false element , then this means that the matrices are not equal to each other.

Since Array2D does not contain methods for filtering or searching, we will implement this ourselves. To do this, we decompose the matrix into a linear list

|> Seq.cast<bool>

And find at least one mismatch

|> Seq.contains false

The Seq.contains function will return true if at least one false value is found in the expanded list . Therefore, we need to invert the result so that the == operator works the way we want

else
    not (matrix1.values
           |> Array2D.mapi (fun x y v -> if matrix2.values.[x, y] <> v then false else true)
           |> Seq.cast<bool>
           |> Seq.contains false)

A matrix O is called a zero or zero matrix if all its elements are equal to zero.

static member O rows cols =
    let array2d = Array2D.zeroCreate rows cols
    { values = array2d }

An example of using this function

let AO = Matrix.O 5 5

I believe that there is nothing complicated that requires explanation, so we continue.

A matrix whose number of rows is equal to the number of columns and equal to n is called a square matrix of order n

Thus, the square matrix has the form.

$$

$$A=\begin{bmatrix}a11&a12&...&a1n\\a21&a22&...&a2n\\...&...&...&...\\an1&an2&...&ann\end{bmatrix}$$

As part of this rule, we will create a function that transforms a rectangular matrix into a square matrix by cutting off all the elements that do not fall into the square.

static member toSquare matrix =

    //    
    let dim = Matrix.sizes matrix

    //   
    let colCount: int = snd dim
    //   
    let rowCount: int = fst dim

    //    
    let length = System.Math.Min (colCount, rowCount)

    //      
    //     
    let zero = Array2D.zeroCreate length length

    //     
    let square = zero |> Array2D.mapi (fun x y _ -> matrix.values.[x, y])

    //   
    { values = square }

Comments in the source code explain how the function works, so let's continue.

A square matrix is ​​called triangular if all its elements below the main diagonal are equal to zero, i.e. the triangular matrix has the form

$$

$$T=\begin{pmatrix}a11&a12&...&a1n\\0&a22&...&a2n\\...&...&...&...\\0&0&...&ann\end{pmatrix}$$

Below is the function code that converts the original matrix into a triangular one. But in our function we will work with a rectangular matrix, that is, it may not be square. The reader can easily modify the function code so that it returns a square triangular matrix using the function that we examined earlier.

static member T matrix =
    let values = matrix.values |> Array2D.mapi (fun x y v -> if y < x then 0 else v)
    { values = values }

The Array2D.mapi function converts the original two-dimensional array into a new one using a handler that takes three parameters

x - row number
y - column number
v - cell contents

if y < x then 0 else v

Here we check whether the element is below the main diagonal, and if so, then fill cell 0. Otherwise, the initial value from the input matrix.

The following is an example of using this function.

let a = array2D [[1;2;3]
                 [4;5;6]
                 [7;8;9]
                 [10;11;12]]
let A = Matrix.ofArray2D a
let R = Matrix.triangular A
printfn "origin = \n %A" A.values
printfn "triangular = \n %A" R.values

We get the following result

origin = 
 [[1; 2; 3]
 [4; 5; 6]
 [7; 8; 9]
 [10; 11; 12]]
triangular = 
 [[1; 2; 3]
 [0; 5; 6]
 [0; 0; 9]
 [0; 0; 0]]

A square matrix is ​​called diagonal if all its elements located outside the main diagonal are equal to zero

static member D matrix =
    let diagonal = matrix.values |> Array2D.mapi (fun x y v -> if x <> y then 0 else v)
    { values = diagonal }

This function is very similar to the previous one, only the verification condition is different. Below is an example of use

let a = array2D [[1;2;3]
                 [4;5;6]
                 [7;8;9]
                 [10;11;12]]
let A = Matrix.ofArray2D a
let R = Matrix.D A
printfn "origin = \n %A" A.values
printfn "diagonal = \n %A" R.values

origin = 
 [[1; 2; 3]
 [4; 5; 6]
 [7; 8; 9]
 [10; 11; 12]]
diagonal = 
 [[1; 0; 0]
 [0; 5; 0]
 [0; 0; 9]
 [0; 0; 0]]

The diagonal matrix is ​​unit and denoted by E if all its elements located on the main diagonal are equal to unity

$$

$$E=\begin{pmatrix}1&0&...&0\\0&1&...&0\\...&...&...&...\\0&0&...&1\end{pmatrix}$$

The implementation of such a matrix on F # looks like this

static member E rows cols =
    let array2d = Array2D.init rows cols (fun x y -> if x = y then 1 else 0)
    { values = array2d }

Matrix Operations with F #

A number of actions can be performed on matrices, as well as on numbers, some of which are similar to operations on numbers, and some are specific.

The sum of two matrices Amn = (aij) and Bmn = (bij) of the same size is the matrix of the same size A + B = Cmn = (cij) , whose elements are equal to the sum of the elements of the matrices A and B located at the corresponding places

$$

$$cij=aij+bij, i=1,2,...,m, j=1,2,...,n$$

Example, for given matrices A and B, we find the sum A + B

$$

$$A=\begin{pmatrix}2&3\\1&-5\\0&6\end{pmatrix}, B=\begin{pmatrix}-3&3\\1&7\\2&0\end{pmatrix}, A+B=\begin{pmatrix}-1&6\\2&2\\2&6\end{pmatrix}$$

Consider the code for adding two matrices

static member (+) (matrix1, matrix2) =
    if Matrix.isEquallySized matrix1 matrix2 then
        let array2d = matrix1.values |> Array2D.mapi (fun x y v -> matrix2.values.[x, y] + v)
        { values = array2d }
    else failwith "matrix1 is not equal to matrix2"

Before adding matrices, you need to make sure that their dimension coincides, otherwise the function throws an exception. Below is an example of using this function

let a = array2D [[2;3]
                 [1;-5]
                 [0;6]]
let A = Matrix.ofArray2D a

let b = array2D [[-3;3]
                 [1;7]
                 [2;0]]
let B = Matrix.ofArray2D b

let R = A+B
printfn "A+B =\n %A" R.values

A+B =
 [[-1; 6]
 [2; 2]
 [2; 6]]

The product of the matrix A = (aij) by the number k is the matrix kA = (kaij) of the same size as the matrix A obtained by multiplying all elements of the matrix A by the number k

Example, for a given matrix A we find matrix 3A

$$

$$A=\begin{pmatrix}2&3&0\\1&5&4\end{pmatrix}, 3A=\begin{pmatrix}6&9&0\\3&15&12\end{pmatrix}$$

static member (*) (value, matrix) = 
    let array2d = matrix.values |> Array2D.mapi (fun _ _ v -> v * value)
    { values = array2d }

Matrix -A = (- 1) * A will be called the opposite matrix A . From this definition, we proceed smoothly to the next

The difference of matrices A and B of equal sizes is the sum of the matrix A and the matrix opposite to B

static member (-) (matrix1: Matrix, matrix2: Matrix) = 
    if Matrix.isEquallySized matrix1 matrix2 then
        matrix1 + (-1)*matrix2
    else failwith "matrix1 is not equal to matrix2"

Two matrices are called consistent if the number of columns of the first is equal to the number of rows of the second

$$

$$A=\begin{pmatrix}2&1\\0&3\end{pmatrix}, B=\begin{pmatrix}0&5&2&1\\1&4&3&7\end{pmatrix}$$

static member isMatched matrix1 matrix2 = 
    let row1Count = matrix1.values.[0,*].Length
    let col2Count = matrix2.values.[*,0].Length

    row1Count = col2Count

Matrix consistency checking is required to multiply them.

The product AB coherent matrix Amn = (aij) and Bnp = (bjk) is the matrix Cmn = (cik) , element cik is calculated as the sum of elements i -th row of the matrix A and corresponding elements of the k th column of the matrix B

Calculate the product of matrices

$$

$$A=\begin{pmatrix}1&0&2\\3&1&0\end{pmatrix}, B=\begin{pmatrix}-1&0\\5&1\\-2&0\end{pmatrix}$$

Decision to determine the product of matrices

$$

$$AB=\begin{pmatrix}1(-1)+0*5+2(-2)&1*0+0*1+2*0\\3(-1)+1*5+0(-2)&3*0+1*1+0*0\end{pmatrix}= \begin{pmatrix}-5&0\\2&1\end{pmatrix}$$

Consider the code for multiplying two matrices

static member (*) (matrix1, (matrix2: Matrix)) =
    if Matrix.isMatched matrix1 matrix2 then
        let row1Count = matrix1.values.[*,0].Length
        let col2Count = matrix2.values.[0,*].Length

        let values = Array2D.zeroCreate row1Count col2Count

        for r in 0..row1Count-1 do
            for c in 0..col2Count-1 do
                let row = Array.toList matrix1.values.[r,*]
                let col = Array.toList matrix2.values.[*,c]

                let cell = List.fold2 (fun acc val1 val2 -> acc + (val1 * val2)) 0 row col
                values.[r,c] <- cell

        { values = values }

    else failwith "matrix1 is not matched to matrix2"

Let's look at the code in more detail.
Before multiplication, make sure that the matrices are consistent

if Matrix.isMatched matrix1 matrix2 then

The resulting matrix will have a dimension in which the number of rows is the same as the first matrix and the number of columns is the same as the second matrix

let row1Count = matrix1.values.[*,0].Length
let col2Count = matrix2.values.[0,*].Length

//        
let values = Array2D.zeroCreate row1Count col2Count

After that, we loop through all the rows and all the columns of the original matrices

for r in 0..row1Count-1 do
    for c in 0..col2Count-1 do
        let row = Array.toList matrix1.values.[r,*]
        let col = Array.toList matrix2.values.[*,c]

We calculate the total value of each cell

let cell = List.fold2 (fun acc val1 val2 -> acc + (val1 * val2)) 0 row col

The List.fold2 function receives two lists at the input (row and column) and passes the following parameters to the handler

acc - the accumulator containing the result of the previous calculation
val1 - the contents of the cell from the first array. In our case, this is the row from the first matrix
val2 - the contents of the cell from the second array, that is, the columns of the second matrix.

Since the matrices are consistent, we are sure that we will not go beyond the arrays.

The handler adds to the accumulator a product of cells from rows and a column and the resulting value will be passed to the next iteration. Thus, the end result of List.fold2 functionwill be the final value of the products of two matrices. It remains only to fill them with the previously created empty matrix

values.[r,c] <- cell

Which will return as a result

{ values = values }

Below is an example of using this function

let a = array2D [[1;0;2]
                 [3;1;0]]
let A = Matrix.ofArray2D a

let b = array2D [[-1;0]
                 [5;1]
                 [-2;0]]
let B = Matrix.ofArray2D b

let R = A*B

printfn "A*B =\n %A" R.values

A1*B1 =
 [[-5; 0]
 [2; 1]]

If k ∈ N, then the kth power of a square matrix A is the product of k matrices A

$$

$$A^k=A*A*...A (k- )$$

Consider the code on F # for the product of a matrix to a power. Tail recursion will be used here in order not to overflow the stack with large powers. Tail recursion is a recursion that the compiler eventually converts to a loop. If possible, it is recommended that you always use tail recursion instead of the usual one, but for this you need each recursion frame to return the final calculated value. This value is usually called the accumulator and is passed to the next recursion frame. That is, unlike regular recursion, which returns the calculated value up the stack, tail recursion passes the calculated value down the stack. Each new recursion frame makes its calculations and adds them to the previously calculated value, which is stored in the battery. After that,As the last recursion frame worked, the accumulator already has a calculated value, which simply returns as a result.

Thus, the compiler can optimize the code and convert it into a regular loop.

//    
// matrix -  
// value -  
static member (^^) (matrix, value) =

    //  ,    
    // m - 
    // p =  
    let inRecPow m p =

        //  
        // acc -  .   Matrix
        // p -     
        //         
        let rec recPow acc p =
            //   
            match p with
            | x when x > 0 ->
                //    
                //      ,      
                let nextAcc = acc*m
                //           
                recPow nextAcc (x-1)
            //    ,    
            | _ -> acc

        //   ,         
        let dim = Matrix.sizes matrix
        let colCount = snd dim
        let rowCount = fst dim

        let u = Matrix.E rowCount colCount
        //           
        recPow u p

    //  ,      
    let powMatrix = inRecPow matrix value
    //   
    { values = powMatrix.values }

The code contains detailed comments on how it works. Requires a little explanation, why is the unit matrix used? It is needed for the first recursion frame and serves as the base value of the accumulator, in which the final result will be accumulated.
Below is an example of using our function

$$

$$A=\begin{pmatrix}1&0\\2&-1\end{pmatrix}$$

$$

$$A^2=AA=\begin{pmatrix}1&0\\2&-1\end{pmatrix}\begin{pmatrix}1&0\\2&-1\end{pmatrix}=\begin{pmatrix}1&0\\0&1\end{pmatrix}$$

$$

$$A^3=AA^2=\begin{pmatrix}1&0\\2&-1\end{pmatrix}\begin{pmatrix}1&0\\0&1\end{pmatrix}=\begin{pmatrix}1&0\\2&-1\end{pmatrix}$$

We calculate the following product

$$

$$f(A)=2A^3-4A^2+3E$$

Where E is the identity matrix. Since we cannot add a number to the matrix, we must add 3E .

$$

$$2\begin{pmatrix}1&0\\2&-1\end{pmatrix}- 4\begin{pmatrix}1&0\\0&1\end{pmatrix}+ 3\begin{pmatrix}1&0\\0&1\end{pmatrix}= \begin{pmatrix}1&0\\4&-3\end{pmatrix}$$

//     
static member (+) (matrix, (value: int)) =
    let dim = Matrix.sizes matrix
    let r = fst dim
    let c = snd dim

    //   
    let unit = Matrix.E r c
    //          
    value*unit + matrix

let a = array2D [[1;0]
                 [2;-1]]
let A = Matrix.ofArray2D a

let R = 2*(A^^3) - 4*(A^^2) + 3
printfn "2*A^3 - 4*A^2 + 3 =\n %A" R.values

2*A5^3 - 4*A5^2 + 3 =
 [[1; 0]
 [4; -3]]

A matrix A ^T whose columns are composed of rows of matrix A with the same numbers and the same sequence of elements is called transposed to matrix A

static member transpose matrix =
    let dim = Matrix.sizes matrix
    let rows = fst dim
    let cols = snd dim

    //      
    let tMatrix = Matrix.O cols rows
    //       
    matrix.values |> Array2D.iteri(fun x y v -> tMatrix.values.[y, x] <- v)

    //  
    tMatrix

Usage example

let a = array2D [[1;2;3]
                 [4;5;6]
                 [7;8;9]
                 [10;11;12]]
let A = Matrix.ofArray2D a
let R6 = Matrix.T A
printfn "origin = \n %A" A.values
printfn "transpose = \n %A" R.values

origin = 
 [[1; 2; 3]
 [4; 5; 6]
 [7; 8; 9]
 [10; 11; 12]]
transpose = 
 [[1; 4; 7; 10]
 [2; 5; 8; 11]
 [3; 6; 9; 12]]

Summary

In this article, we examined examples of the implementation and use of matrices from the theory of linear algebra. As well as the basic mathematical operations on them, using a functional approach in F #. I hope the reader can appreciate the flexibility that functional languages ​​provide.

The full source code of the matrix module, fragments of which were considered in the framework of the article, you can find on the github.

Github Matrix.fs