ParaMonte Fortran 2.0.0 Parallel Monte Carlo and Machine Learning LibrarySee the latest version documentation.
pm_matrixInv Module Reference

This module contains abstract and concrete derived types and procedures related to the inversion of square matrices.
More...

## Data Types

interface  getMatInv
Generate and return the full inverse of an input matrix of general or triangular form directly or through its input the LU/Cholesky factorization.
More...

type  inversion_type
This is a concrete derived type whose instances are exclusively used to request inversion operation on a given matrix within an interface of a procedure of the ParaMonte library.
More...

interface  setMatInv
Generate and return the full inverse of a general or triangular matrix or a subset of the inverse of a positive-definite matrix complementary to its specified Cholesky factorization subset.
More...

## Variables

character(*, SK), parameter MODULE_NAME = "@pm_matrixInv"

type(inversion_type), parameter inversion = inversion_type()
This is a scalar parameter object of type inversion_type that is exclusively used to request no transpose of a given matrix within an interface of a procedure of the ParaMonte library.
More...

## Detailed Description

This module contains abstract and concrete derived types and procedures related to the inversion of square matrices.

# Inversion operation

In linear algebra, an $$n$$-by- $$n$$ square matrix $$A$$ is called invertible (also nonsingular, nondegenerate), if there exists an $$n$$-by- $$n$$ square matrix $$B$$ such that,

$$\mathbf{AB} = \mathbf{BA} = \mathbf{I}_{n}~,$$

where $$I_n$$ denotes the $$n$$-by- $$n$$ identity matrix and the multiplication used is ordinary matrix multiplication.
If this is the case, then the matrix $$B$$ is uniquely determined by $$A$$, and is called the **(multiplicative) inverse** of $$A$$, denoted by $$A^{−1}$$.
source inversion is the process of finding the matrix $$B$$ that satisfies the prior equation for a given invertible matrix $$A$$.
Over a field, a square matrix that is not invertible is called singular or degenerate.
A square matrix with entries in a field is singular if and only if its determinant is zero.
Singular matrices are rare in the sense that if a square matrix entries are randomly selected from any bounded region on the number line or complex plane, the probability that the matrix is singular is $$0$$, that is, it will almost never be singular.
Non-square matrices do not have an inverse.
However, in some cases such a matrix may have a left inverse or right inverse.
If $$A$$ is $$m$$-by- $$n$$ and the rank of $$A$$ is equal to $$n$$ ( $$n \leq m$$), then $$A$$ has a left inverse, an $$n$$-by- $$m$$ matrix $$B$$ such that $$BA = I_n$$.
If $$A$$ has rank $$m$$ ( $$m \leq n$$), then it has a right inverse, an $$n$$-by- $$m$$ matrix $$B$$ such that $$AB = I_m$$.

## Inverse matrix properties

The following properties hold for an invertible matrix $$A$$:

1. $$(\mathbf{A}^{-1})^{-1} = \mathbf{A}$$.
2. $$(k\mathbf{A})^{-1} = k^{-1}\mathbf{A}^{-1}$$ for nonzero scalar $$k$$.
3. $$(\mathbf{Ax})^{+} = \mathbf{x}^{+}\mathbf{A}^{-1}$$ if $$A$$ has orthonormal columns, where $$+$$ denotes the Moore–Penrose inverse and $$x$$ is a vector.
4. $$(\mathbf{A}^{\mathrm{T}})^{-1} = (\mathbf{A}^{-1})^{\mathrm{T}}$$.
5. For any invertible $$n$$-by- $$n$$ matrices $$A$$ and $$B$$, $$(\mathbf{AB})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}$$.
More generally, if $$\mathbf{A}_1, \dots, \mathbf{A}_{k}$$ are invertible $$n$$-by- $$n$$ matrices, then $$(\mathbf{A}_{1}\mathbf{A}_{2} \cdots \mathbf{A}_{k-1}\mathbf{A}_{k})^{-1} = \mathbf{A}_{k}^{-1}\mathbf{A}_{k-1}^{-1} \cdots \mathbf{A}_{2}^{-1}\mathbf{A}_{1}^{-1}$$.
6. $$\det\mathbf{A}^{-1} = (\det \mathbf{A})^{-1}$$.
7. The rows of the inverse matrix $$V$$ of a matrix $$U$$ are orthonormal to the columns of $$U$$.

## Inverse matrix computation

The common approach to computing the inverse matrix stems from its definition.
The inverse matrix $$A^{-1}$$ of a square matrix $$A$$ is a square matrix such that $$AA^{-1} = I$$, where $$I$$ is the identity matrix.
Depending on the class of the square matrix $$A$$, there are several approaches that can be taken to compute its inverse.
However, all such methods attempt to factorize the matrix first (for example, Cholesky decomposition or LU decomposition) and cast the problem into seeking the solution to a system of linear equations.
For example, for a general square matrix of shape $$3\times 3$$, the corresponding system of equations to solve would be:

$$\begin{bmatrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{bmatrix} \begin{bmatrix} A_{11}^{-1} & A_{12}^{-1} & A_{13}^{-1} \\ A_{21}^{-1} & A_{22}^{-1} & A_{23}^{-1} \\ A_{31}^{-1} & A_{32}^{-1} & A_{33}^{-1} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} ~,$$

where the second matrix on the left hand side is the inverse of the first square matrix $$A$$.
The inverse matrix can be constructed as the collection of solutions to $$n = 3$$ systems of equations of the form $$Ax=b$$ with different right hand sides $$b$$ matrices and $$x$$ representing $$n^{\mathrm{th}}$$ column of the inverse matrix.
The first system of equations to solve for the above problem would be:

$$\begin{bmatrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{bmatrix} \begin{bmatrix} A_{11}^{-1} \\ A_{21}^{-1} \\ A_{31}^{-1} \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \\ \end{bmatrix} ~.$$

The second column of the inverse can be computed by changing $$b$$ to $$[0,1,0]^T$$, the third column with $$[0,0,1]^T$$, and so on.
The task of computing the inverse of the matrix is now reduced to solving a series of systems of linear equations.
This can be done if the matrix $$A$$ is factorized into lower and upper triangular matrices,

$$\begin{bmatrix} A_{11} & A_{12} & A_{13} \\ A_{21} & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{bmatrix} = \begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} ~, \end{bmatrix}$$

such that the above system can be written as,

$$\begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \end{bmatrix} \left( \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \end{bmatrix} \begin{bmatrix} A_{11}^{-1} \\ A_{21}^{-1} \\ A_{31}^{-1} \end{bmatrix} \right) = \begin{bmatrix} 1 \\ 0 \\ 0 ~. \end{bmatrix}$$

The above form of equations implied that the two systems of triangular equations have to be solved to obtain one column of the inverse matrix.
However, this method is fast because only back- and forward-substitution are required to solve for the column vectors after the initial factorization of the matrix $$A$$.
The most common factorization methods used for $$A$$ to obtain its inverse are:

1. Pivoted LU factorization for general square matrices.
2. Cholesky factorization for Symmetric/Hermitian Positive-Definite square matrices.

# Inversion-Transposition operation

There are also special matrix operations that mix inversion with Symmetric and Hermitian each having corresponding matrix classes:**

1. Orthogonal transposition
A square matrix whose transpose is equal to its inverse is called an Orthogonal matrix.
In other words, $$A$$ is Orthogonal if $$\mathbf{A}^{\up{T}} = \mathbf{A}^{-1}$$.
The corresponding transposition is called Orthogonal denoted by the operator $$\cdot^{\up{-T}}$$.
2. Unitary transposition
A square complex matrix whose transpose is equal to its conjugate inverse is called a Unitary matrix.
In other words, $$A$$ is Unitary if $$\mathbf{A}^{\up{T}} = {\overline{\mathbf{A}^{-1}}}$$.
The corresponding transposition is called Unitary denoted by the operator $$\cdot^{\up{-H}}$$.
pm_matrixDet
pm_matrixLUP
pm_matrixChol
pm_matrixTrans
Benchmarks:

Benchmark :: The runtime performance of setMatInv for various methods and inverse matrix subsets.

1! Test the performance of Cholesky factorization computation using an assumed-shape interface vs. explicit-shape interface.
2program benchmark
3
4 use pm_kind, only: IK, LK, RKG => RKD, SK
5 use pm_matrixCopy, only: setMatCopy, rdpack, uppDia, lowDia, transHerm
6 use pm_distUnif, only: rngx_type => xoshiro256ssw_type
8 use pm_distUnif, only: getUnifRand
9 use pm_bench, only: bench_type
10
11 implicit none
12
13 integer(IK) :: itry, ntry
14 integer(IK) :: i
15 integer(IK) :: iarr
16 integer(IK) :: fileUnit
17 integer(IK) , parameter :: NARR = 10_IK
18 integer(IK) , allocatable :: rperm(:)
19 real(RKG) , allocatable :: mat(:,:), inv(:,:)
20 type(bench_type), allocatable :: bench(:)
21 integer(IK) , parameter :: nsim = 2**NARR
22 integer(IK) :: rank
23 real(RKG) :: dumm
24 type(rngx_type) :: rngx
25
26 rngx = rngx_type()
27 bench = [ bench_type(name = SK_"setMatInvLow", exec = setMatInvLow, overhead = setOverhead) &
28 , bench_type(name = SK_"setMatInvUpp", exec = setMatInvUpp, overhead = setOverhead) &
29 , bench_type(name = SK_"setMatInvLUP", exec = setMatInvLUP, overhead = setOverhead) &
30#if LAPACK_ENABLED
31 , bench_type(name = SK_"lapack_dpotri", exec = lapack_dpotri, overhead = setOverhead) &
32#endif
33 ]
34
35 write(*,"(*(g0,:,' '))")
36 write(*,"(*(g0,:,' '))") "inverse matrix benchmarking..."
37 write(*,"(*(g0,:,' '))")
38
39 open(newunit = fileUnit, file = "main.out", status = "replace")
40
41 write(fileUnit, "(*(g0,:,','))") "rank", (bench(i)%name, i = 1, size(bench))
42
43 dumm = 0._RKG
44 loopOverMatrixSize: do iarr = 1, NARR
45
46 rank = 2**iarr
47 ntry = nsim / rank
48 call setResized(rperm, rank)
49 call setResized(inv, [rank, rank])
50 mat = getUnifRand(1._RKG, 2._RKG, rank, rank + 1_IK)
51 write(*,"(*(g0,:,' '))") "Benchmarking setMatInv() algorithms with array size", rank, ntry
52
53 do i = 1, size(bench)
54 bench(i)%timing = bench(i)%getTiming()
55 end do
56
57 write(fileUnit,"(*(g0,:,','))") rank, (bench(i)%timing%mean / ntry, i = 1, size(bench))
58
59 end do loopOverMatrixSize
60 write(*,"(*(g0,:,' '))") dumm
61
62 close(fileUnit)
63
64contains
65
66 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 ! procedure wrappers.
68 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
69
71 do itry = 1, ntry
72 call setMatrix()
73 end do
74 end subroutine
75
76 subroutine setMatrix()
77 !integer(IK) :: i
78 !call random_number(mat)
79 !mat = mat * 1.e-4_RKG
80 !do i = 1, size(mat, dim = 1, kind = IK)
81 ! mat(i,i+1) = 1._RKG
82 !end do
83 !use pm_distCov, only: setCovRand
84 !call setCovRand(rngx, mat(:,2:rank+1)) ! causes numeric overflow for large matrix ranks.
85 use pm_matrixInit, only: setMatInit, uppLowDia
86 call setMatInit(mat(:,2:rank+1), uppLowDia, 0._RKG, 0._RKG, 1._RKG)
87 dumm = dumm - mat(rank, rank) - mat(rank-1, rank-1)
88 end subroutine
89
90#if LAPACK_ENABLED
91 subroutine lapack_dpotri()
92 integer(IK) :: info
93 do itry = 1, ntry
94 call setMatrix()
95 call dpotrf("U", rank, mat(:,2:rank+1), rank, info)
96 if (info /= 0_IK) error stop
97 call dpotri("U", rank, mat(:,2:rank+1), rank, info)
98 if (info /= 0_IK) error stop
99 call setMatCopy(mat(:,2:rank+1), rdpack, mat(:,1:rank), rdpack, uppDia, transHerm) ! symmetrize inverse matrix.
100 !dumm = dumm + mat(rank,rank) + mat(rank-1,rank-1)
101 end do
102 end subroutine
103#endif
104
105 subroutine setMatInvLow()
106 use pm_matrixInv, only: setMatInv, choUpp
107 use pm_matrixChol, only: setMatChol, nothing
108 integer(IK) :: info
109 do itry = 1, ntry
110 call setMatrix()
111 call setMatChol(mat(:,2:rank+1), uppDia, info, mat(:,2:rank+1), nothing)
112 if (info /= 0_IK) error stop
113 call setMatInv(mat(:,1:rank), mat(:,2:rank+1), choUpp)
114 !dumm = dumm + mat(rank,rank) + mat(rank-1,rank-1)
115 end do
116 end subroutine
117
118 subroutine setMatInvUpp()
119 use pm_matrixInv, only: setMatInv, choLow, choUpp
120 use pm_matrixChol, only: setMatChol
121 integer(IK) :: info
122 do itry = 1, ntry
123 call setMatrix()
124 call setMatChol(mat(:,2:rank+1), uppDia, info, mat(:,1:rank), transHerm)
125 if (info /= 0_IK) error stop
126 call setMatInv(mat(:,2:rank+1), mat(:,1:rank), choLow)
127 !dumm = dumm + mat(rank,rank) + mat(rank-1,rank-1)
128 end do
129 end subroutine
130
131 subroutine setMatInvLUP()
132 use pm_matrixLUP, only: setMatLUP
133 use pm_matrixInv, only: setMatInv
134 integer(IK), parameter :: offset = 1
135 integer(IK) :: info
136 do itry = 1, ntry
137 call setMatrix()
138 call setMatLUP(mat(:,2:rank+1), rperm, info)
139 if (info /= 0_IK) error stop
140 call setMatInv(inv, mat(:,2:rank+1), rperm)
141 !dumm = dumm + mat(rank,rank) + mat(rank-1,rank-1)
142 end do
143 end subroutine
144
145end program benchmark
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Generate and return an object of type timing_type containing the benchmark timing information and sta...
Definition: pm_bench.F90:574
Generate and return a scalar or a contiguous array of rank 1 of length s1 of randomly uniformly distr...
Compute and return the lower/upper-triangle of the Cholesky factorization of the input Symmetric/Her...
Copy a desired subset of the input source matrix of arbitrary shape (:) or (:,:) to the target subset...
Set the upper/lower triangle and the diagonal elements of the input matrix of arbitrary shape (:,...
Generate and return the full inverse of a general or triangular matrix or a subset of the inverse of ...
Return the LU-Pivoted decomposition of the input square matrix mat(ndim,ndim).
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
This module contains abstract interfaces and types that facilitate benchmarking of different procedur...
Definition: pm_bench.F90:41
This module contains classes and procedures for computing various statistical quantities related to t...
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter RKD
The double precision real kind in Fortran mode. On most platforms, this is an 64-bit real kind.
Definition: pm_kind.F90:568
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
This module contains procedures and generic interfaces for computing the Cholesky factorization of po...
This module contains procedures and generic interfaces relevant to copying (diagonal or upper/lower t...
This module contains procedures and generic interfaces relevant to generating and initializing matric...
This module contains abstract and concrete derived types and procedures related to the inversion of s...
This module contains procedures and generic interfaces relevant to the partially LU Pivoted decomposi...
This is the class for creating benchmark and performance-profiling objects.
Definition: pm_bench.F90:386
This is the derived type for declaring and generating objects of type xoshiro256ssw_type containing a...
subroutine bench(sort, arraySize)

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Postprocessing of the benchmark output
1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4import pandas as pd
5import numpy as np
6
7import os
8dirname = os.path.basename(os.getcwd())
9
10fontsize = 14
11
12df = pd.read_csv("main.out", delimiter = ",")
13colnames = list(df.columns.values)
14
15
18
19ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
20ax = plt.subplot()
21
22for colname in colnames[1:]:
23 plt.plot( df[colnames[0]].values
24 , df[colname].values
25 , linewidth = 2
26 )
27
28plt.xticks(fontsize = fontsize)
29plt.yticks(fontsize = fontsize)
30ax.set_xlabel(colnames[0], fontsize = fontsize)
31ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
32ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
33ax.set_xscale("log")
34ax.set_yscale("log")
35plt.minorticks_on()
36plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
37ax.tick_params(axis = "y", which = "minor")
38ax.tick_params(axis = "x", which = "minor")
39ax.legend ( colnames[1:]
40 #, loc='center left'
41 #, bbox_to_anchor=(1, 0.5)
42 , fontsize = fontsize
43 )
44
45plt.tight_layout()
46plt.savefig("benchmark." + dirname + ".runtime.png")
47
48
51
52ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
53ax = plt.subplot()
54
55plt.plot( df[colnames[0]].values
56 , np.ones(len(df[colnames[0]].values))
57 , linestyle = "--"
58 #, color = "black"
59 , linewidth = 2
60 )
61for colname in colnames[2:]:
62 plt.plot( df[colnames[0]].values
63 , df[colname].values / df[colnames[1]].values
64 , linewidth = 2
65 )
66
67plt.xticks(fontsize = fontsize)
68plt.yticks(fontsize = fontsize)
69ax.set_xlabel(colnames[0], fontsize = fontsize)
70ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
71ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}().".format(colnames[1]), fontsize = fontsize)
72ax.set_xscale("log")
73ax.set_yscale("log")
74plt.minorticks_on()
75plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
76ax.tick_params(axis = "y", which = "minor")
77ax.tick_params(axis = "x", which = "minor")
78ax.legend ( colnames[1:]
79 #, bbox_to_anchor = (1, 0.5)
80 #, loc = "center left"
81 , fontsize = fontsize
82 )
83
84plt.tight_layout()
85plt.savefig("benchmark." + dirname + ".runtime.ratio.png")

Visualization of the benchmark output

Benchmark moral
1. The procedures under the generic interface setMatInv use an unblocked approach to computing the matrix inverse.
However, specifying an upper-triangular Cholesky factor along with lower-triangle for the matrix inverse can potentially result in faster calculations as all matrix operations within the algorithm become column-major meaning that all memory access become local.
Test:
test_pm_matrixInv

Final Remarks

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.

1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin

## ◆ inversion

 type(inversion_type), parameter pm_matrixInv::inversion = inversion_type()

This is a scalar parameter object of type inversion_type that is exclusively used to request no transpose of a given matrix within an interface of a procedure of the ParaMonte library.

For example usage, see the documentation of the target procedure requiring this object.

nothing
trans

Final Remarks

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.

1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin

Definition at line 277 of file pm_matrixInv.F90.

## ◆ MODULE_NAME

 character(*,SK), parameter pm_matrixInv::MODULE_NAME = "@pm_matrixInv"

Definition at line 231 of file pm_matrixInv.F90.