Source Codes
MATRIX INVERSION
Private Sub Matrix_Inversion(ByVal A(,) As Single, ByVal Bound As Integer, ByRef A_inv(,) As Single) '---- Gauss Jordan Inverse Matrix |A I| = |I B| where B = A^-1 Dim _triang(Bound, (Bound * 2) + 1) As Single Dim p, q As Integer '----Forming |A I| matrix For p = 0 To Bound For q = 0 To Bound _triang(p, q) = A(p, q) If p = q Then _triang(p, q + (Bound + 1)) = 1 End If Next Next '----Gauss Jordan for reducing to |B I| Dim line_1 As Integer Dim temp, pivot As Single For k = 0 To Bound 'Bring a non-zero element first by changes lines if necessary If _triang(k, k) = 0 Then For n = k To Bound If _triang(n, k) <> 0 Then line_1 = n : Exit For 'Finds line_1 with non-zero element Next n 'Change line k with line_1 For m = k To ((Bound * 2) + 1) temp = _triang(k, m) _triang(k, m) = _triang(line_1, m) _triang(line_1, m) = temp Next m End If pivot = _triang(k, k) For n = k To (Bound * 2) + 1 _triang(k, n) = _triang(k, n) / pivot Next n 'For other lines, make a zero element by using: 'Ai1=Aij-A11*(Aij/A11) 'and change all the line using the same formula for other elements Dim mult_1 As Single For n = 0 To Bound If n = k And n = Bound Then Exit For If n = k And n < Bound Then n = n + 1 'Do not change that element (already equals to 1), go for next one If _triang(n, k) <> 0 Then 'if it is zero, stays as it is mult_1 = _triang(n, k) / _triang(k, k) For m = k To (Bound * 2) + 1 _triang(n, m) = _triang(n, m) - (_triang(k, m) * mult_1) Next m End If Next n Next k For p = 0 To Bound For q = 0 To Bound A_inv(p, q) = _triang(p, q + (Bound + 1)) Next Next End Sub Private Sub Matrix_Transpose(ByVal A(,) As Single, ByVal Bound As Integer, ByRef AT(,) As Single) Dim i, j As Integer For i = 0 To Bound For j = 0 To Bound AT(i, j) = A(j, i) Next Next End SubCHOLESKI DECOMPOSITION
CHOLESKI DECOMPOSITION
Private Sub Choleski_Decomposition(ByVal A(,) As Single, ByVal n As Integer, ByRef L(,) As Single) Dim Sum1, Sum2, Sum3 As Single ReDim L(n, n) L(0, 0) = Math.Sqrt(A(0, 0)) For j = 1 To n L(j, 0) = A(j, 0) / L(0, 0) Next For i = 1 To n - 1 For k = 0 To i Sum1 = Sum1 + Math.Pow(L(i, k), 2) Next L(i, i) = Math.Sqrt(A(i, i) - Sum1) For j = i + 1 To n For k = 0 To i Sum2 = Sum2 + (L(j, k) * L(i, k)) Next L(j, i) = ((A(j, i) - Sum2) / (L(i, i))) Next Next For k = 0 To n - 1 Sum3 = Sum3 + Math.Pow(L(n, k), 2) Next L(n, n) = Math.Sqrt(A(n, n) - Sum3) End SubGAUSS ELIMINATION
GAUSS ELIMINATION
#Region "Gauss Elimination Method" '-----Redefined Gauss Elimination Procedure Private Sub Gauss(ByVal A(,) As Single, ByVal B() As Single, ByRef X() As Single, ByVal Bound As Integer) '___________________________________________________________________________________________ '---------------- HouseHolder Transformation developed by Samson Mano ---------------------- '-------------------- based on the following notes ----------------------------------------- '----------------- 1. https://en.wikipedia.org/wiki/Gaussian_elimination-------------------- '___________________________________________________________________________________________ Dim Triangular_A(Bound, Bound + 1) As Single Dim soln(Bound) As Single 'Solution matrix For m = 0 To Bound For n = 0 To Bound Triangular_A(m, n) = A(m, n) Next Next '.... substituting the force to triangularmatrics.... For n = 0 To Bound Triangular_A(n, Bound + 1) = B(n) Next ForwardSubstitution(Triangular_A, Bound) ReverseElimination(Triangular_A, X, Bound) End Sub Private Sub ForwardSubstitution(ByRef _triang(,) As Single, ByVal bound As Integer) 'Forward Elimination 'Dim _fraction As Single For k = 0 To bound - 1 For i = k + 1 To bound If _triang(k, k) = 0 Then Continue For End If _triang(i, k) = _triang(i, k) / _triang(k, k) For j = k + 1 To bound + 1 _triang(i, j) = _triang(i, j) - (_triang(i, k) * _triang(k, j)) Next Next Next End Sub Private Sub ReverseElimination(ByRef _triang(,) As Single, ByRef X() As Single, ByVal bound As Integer) 'Back Substitution For i = 0 To bound X(i) = _triang(i, bound + 1) Next For i = bound To 0 Step -1 For j = i + 1 To bound X(i) = X(i) - (_triang(i, j) * X(j)) Next X(i) = X(i) / _triang(i, i) Next End Sub#End RegionEIGEN VALUES & EIGEN VECTORS USING HOUSEHOLDER TRANSFORMATION/ QR DECOMPOSITION
EIGEN VALUES & EIGEN VECTORS USING HOUSEHOLDER TRANSFORMATION/ QR DECOMPOSITION
Householder transformations are widely used in numerical linear algebra, to perform QR decompositions and in the first step of the QR algorithm. The QR algorithm for the computation of eigenvalues, which is based on the QR-decomposition, is considered to be one of the 10 most important algorithms of the 20th century
#Region "EigenValues and EigenVectors Calculator Using HouseHolder Transformation and QR Decomposition" '___________________________________________________________________________________________________________________________________________________________________________________________ '___________________________________________________________________________________________________________________________________________________________________________________________ Private Sub EigenValues_and_EigenVectors_HouseHolder_Transfromation(ByVal A_Matrix(,) As Double, ByVal n As Integer, ByRef EigenValues() As Double, ByRef Eigen_Vector(,) As Double) '___________________________________________________________________________________________ '---------------- EigenValues and EigenVector calculator developed by Samson Mano ---------- '------------ based on HouseHolder Transfromation and QR Decomposition --------------------- '-------------- Eigen values and Eigen Vectors of A_Matrix is calculated ------------------ '--------- for 4x4 matrix n = 3 , 7x7 matrix n =4 and ExE matrix n = E-1 !!!! Important Note '___________________________________________________________________________________________ Dim TPPk_Matrix(n, n) As Double HouseHolder_Transformation(A_Matrix, n, TPPk_Matrix) Dim Possib_EigenVector(n, n) As Double Dim lv1, lv2, lv3 As Integer '----- keep unit matrix for possib_Eigenvector For lv1 = 0 To n Step +1 Possib_EigenVector(lv1, lv1) = 1 Next Do Dim TQMatrix(n, n) As Double Dim TRMatrix(n, n) As Double QRDecomposition(A_Matrix, n, TQMatrix, TRMatrix) '---- Update A matrix to A1 which is A = QR => Q^-1A = R => Q^-1AQ = RQ = A1 and Q^-1 = Q^T For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 A_Matrix(lv1, lv2) = 0 For lv3 = 0 To n Step +1 A_Matrix(lv1, lv2) = A_Matrix(lv1, lv2) + (TRMatrix(lv1, lv3) * TQMatrix(lv3, lv2)) Next Next Next '---- Update EigenVector here which is Q1*Q2*Q3 .... Qk Dim Temp_Q(n, n) As Double For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 Temp_Q(lv1, lv2) = 0 For lv3 = 0 To n Step +1 Temp_Q(lv1, lv2) = Temp_Q(lv1, lv2) + (Possib_EigenVector(lv1, lv3) * TQMatrix(lv3, lv2)) Next Next Next For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 Possib_EigenVector(lv1, lv2) = Temp_Q(lv1, lv2) Next Next '---- Exit if lower triangle of A_Matrix becomes Zero If Check_IfLowerTriangleisZero(A_Matrix, n) = True Then Exit Do End If Loop '-------------- Fixing the EigenValues here For lv1 = 0 To n EigenValues(lv1) = A_Matrix(lv1, lv1) Next '------------ Fixing the EigenVectors here For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 Eigen_Vector(lv1, lv2) = 0 For lv3 = 0 To n Step +1 Eigen_Vector(lv1, lv2) = Eigen_Vector(lv1, lv2) + (TPPk_Matrix(lv1, lv3) * Possib_EigenVector(lv3, lv2)) Next Next Next '********** Sorting EigenValues / EigenVectors For lv1 = 0 To n - 1 Step +1 For lv2 = lv1 + 1 To n Step +1 'Sorting Ascending Minimum to Maximum If EigenValues(lv1) > EigenValues(lv2) Then SortSwap(EigenValues(lv1), EigenValues(lv2)) For lv3 = 0 To n Step +1 SortSwap(Eigen_Vector(lv3, lv1), Eigen_Vector(lv3, lv2)) Next End If Next Next End Sub Private Sub SortSwap(ByRef a As Double, ByRef b As Double) Dim temp As Double temp = a a = b b = temp End Sub Private Function Check_IfLowerTriangleisZero(ByRef A_Matrix(,) As Double, ByVal n As Integer) As Boolean Dim lv1, lv2 As Integer Dim Chk_Bool As Boolean Chk_Bool = True For lv1 = 1 To n Step +1 For lv2 = 0 To (lv1 - 1) Step +1 If Math.Round(A_Matrix(lv1, lv2), 5) <> 0 Then Chk_Bool = False GoTo Ex1 End If Next NextEx1: Return Chk_Bool End Function Private Sub HouseHolder_Transformation(ByRef A_Matrix(,) As Double, ByVal n As Integer, ByRef PPk_Matrix(,) As Double) '___________________________________________________________________________________________ '---------------- HouseHolder Transformation developed by Samson Mano ---------------------- '-------------------- based on the following notes ----------------------------------------- '----------------- 1. https://en.wikipedia.org/wiki/Householder_transformation-------------- '___________________________________________________________________________________________ Dim k As Integer '--- Fix Unit Matrix for P0 so that we can multiply p1p2..pk For k = 0 To n Step +1 PPk_Matrix(k, k) = 1 Next For k = 0 To n - 2 Step +1 '----- Find the alpha value alpha = -1*(sgn(a(k+1,k))*sqrt(a(j,k)^2) ->j from k+1 to n Dim alpha As Double = 0 Dim lv1, lv2, lv3 As Integer For lv1 = k + 1 To n Step +1 alpha = alpha + (A_Matrix(lv1, k) ^ 2) Next alpha = (-1) * (Math.Sign(A_Matrix(k + 1, k))) * Math.Sqrt(alpha) '----- Find the r value r_value = math.sqrt((1/2)*(alpha^2 - (A_matrix(k+1,k)*alpha)) Dim r_value As Double r_value = Math.Sqrt((1 / 2) * ((alpha ^ 2) - (A_Matrix(k + 1, k) * alpha))) '----- Find the v matrix v(0)=v(1)=v(2)=...=v(k)=0 and v(k+1)=(a(k+1,k)-alpha)/(2*r) and v(j)=a(j,k)/(2*r) -> j = k+2,k+3,...,n Dim v_matrix(n) As Double For lv1 = 0 To k Step +1 v_matrix(lv1) = 0 Next v_matrix(k + 1) = (A_Matrix(k + 1, k) - alpha) / (2 * r_value) For lv1 = k + 2 To n Step +1 v_matrix(lv1) = A_Matrix(lv1, k) / (2 * r_value) Next '----- Find the vvT matrix and P_Matrix Dim vvT(n, n) As Double Dim I_Matix(n, n) As Double 'Multiply v matrix with vT matrix For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 vvT(lv1, lv2) = 0 I_Matix(lv1, lv2) = 0 vvT(lv1, lv2) = v_matrix(lv1) * v_matrix(lv2) Next I_Matix(lv1, lv1) = 1 '--- All the diagonals are 1 Next 'Multiply vvT matrix with 2 For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 vvT(lv1, lv2) = 2 * vvT(lv1, lv2) Next Next 'Forming P_Matrix here Dim P_Matrix(n, n) As Double For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 P_Matrix(lv1, lv2) = I_Matix(lv1, lv2) - vvT(lv1, lv2) Next Next '----- Store P_matrix as mutiples to use in transformation PPk = p1p2p3...pk Dim Temp_P(n, n) As Double For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 Temp_P(lv1, lv2) = 0 For lv3 = 0 To n Step +1 Temp_P(lv1, lv2) = Temp_P(lv1, lv2) + (PPk_Matrix(lv1, lv3) * P_Matrix(lv3, lv2)) Next Next Next For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 PPk_Matrix(lv1, lv2) = Temp_P(lv1, lv2) Next Next '----- Find P_Matrix*A_Matrix*P_Matrix which is the A_Matrix HouseHolder Transformation Dim AP_Matrix(n, n) As Double 'Multiply A_Matrix with P_Matrix to find intermittent AP_Matrix For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 AP_Matrix(lv1, lv2) = 0 For lv3 = 0 To n Step +1 AP_Matrix(lv1, lv2) = AP_Matrix(lv1, lv2) + (A_Matrix(lv1, lv3) * P_Matrix(lv3, lv2)) Next Next Next 'Multiply P_Matrix with AP_Matrix For lv1 = 0 To n Step +1 For lv2 = 0 To n Step +1 A_Matrix(lv1, lv2) = 0 For lv3 = 0 To n Step +1 A_Matrix(lv1, lv2) = A_Matrix(lv1, lv2) + (P_Matrix(lv1, lv3) * AP_Matrix(lv3, lv2)) Next Next Next Next End Sub Private Sub QRDecomposition(ByRef A_Matrix(,) As Double, ByRef n As Integer, ByRef Q_Matrix(,) As Double, ByRef R_Matrix(,) As Double) '___________________________________________________________________________________________ '---------------- QR Algorithgm developed by Samson Mano ----------------------------------- '-------------------- based on the following notes ----------------------------------------- '----------------- 1. http://ranger.uta.edu/~weems/NOTES4351/hansen.pdf--------------------- '----------------- 2. http://web.engr.illinois.edu/~heath/scicomp/notes/chap03.pdf---------- '----------------- 3. http://www.math.usm.edu/lambers/mat610/sum10/lecture9.pdf------------- '___________________________________________________________________________________________ '--------------------Initialization A- Matrix Dim A(n, n) As Double Dim R(n - 1)(,) As Double Dim i, j, k As Integer For i = 0 To n For j = 0 To n A(i, j) = A_Matrix(i, j) Next Next '___________________________________________________________________________ For k = 0 To (n - 1) Step 1 Dim A_Norm As Double = 0 For i = k To n A_Norm = A_Norm + (A(i, k) ^ 2) '---- Finding the length of column vector Next A_Norm = Math.Sqrt(A_Norm) Dim Alpha_A As Double = If(A_Norm < 0, A_Norm, -A_Norm) '----- Forming V - Matrix v = a - alpha * (e) Dim V(n) As Double For i = k To n Step 1 V(i) = A(i, k) Next V(k) = V(k) - Alpha_A '------------------------- '--- H = I - (2/V^T*V)* V *V^T V_VT_Mul(V, k, n, R) For i = k To n Step 1 Dim V_Ratio As Double Dim A_Col(n) As Double For j = k To n Step 1 A_Col(j) = A(j, i) '---- Transformation for every column Next V_Ratio = VT_A_Mul(V, A_Col, k, n) * VT_V_Mul(V, k, n) For j = k To n Step 1 A(j, i) = A_Col(j) - (V_Ratio * V(j)) '---- Transformation for every column Next Next Next '------------- Multiplying all the H - Matrix to form Q_Matrix Dim Q(0)(,) As Double ReDim Q(0)(n, n) For i = 0 To n Q(0)(i, i) = 1 '----- Creating I - matrix Next For k = 0 To (n - 1) Step 1 H_Multiplication(Q, R, k, n, Q) Next '_______________________________________________________________ 'Transferring the result For i = 0 To n For j = 0 To n R_Matrix(i, j) = A(i, j) Next Next For i = 0 To n For j = 0 To n Q_Matrix(i, j) = Q(0)(i, j) Next Next End Sub Private Sub H_Multiplication(ByRef A()(,) As Double, _ ByRef B()(,) As Double, _ ByRef k As Integer, _ ByRef n As Integer, _ ByRef C()(,) As Double) '------- Ordinary Matrix multiplication Dim T(n, n) As Double For i = 0 To n For j = 0 To n T(i, j) = 0 For d = 0 To n T(i, j) = T(i, j) + (A(0)(i, d) * B(k)(d, j)) Next Next Next ReDim C(0)(n, n) For i = 0 To n For j = 0 To n C(0)(i, j) = T(i, j) Next Next End Sub Private Sub V_VT_Mul(ByRef V() As Double, ByRef k As Integer, ByRef n As Integer, ByRef H()(,) As Double) ReDim H(k)(n, n) Dim c As Double = VT_V_Mul(V, k, n) '------- Initilaizing For i = 0 To n For j = 0 To n If i = j Then H(k)(i, j) = 1 Else H(k)(i, j) = 0 End If Next Next '---- Doing the operation H = I - (2/vT*V)*(V*VT) For i = k To n For j = k To n If i = j Then H(k)(i, j) = 1 - (c * V(i) * V(j)) Else H(k)(i, j) = 0 - (c * V(i) * V(j)) End If Next Next End Sub Private Function VT_A_Mul(ByRef V() As Double, ByRef A() As Double, ByRef k As Integer, ByRef n As Integer) '-------- Function calculates ratio of |V^T*A/V^T*V| Dim VtV As Double = 0 ' Holds the value of V^T*V Dim VtA As Double = 0 ' Holds the value of V^T*A For i = k To n VtA = VtA + (V(i) * A(i)) Next Return (VtA) End Function Private Function VT_V_Mul(ByRef V() As Double, ByRef k As Integer, ByRef n As Integer) '-------- Function calculates ratio of |V^T*A/V^T*V| Dim VtV As Double = 0 ' Holds the value of V^T*V Dim VtA As Double = 0 ' Holds the value of V^T*A For i = k To n VtV = VtV + (V(i) ^ 2) Next Return (2 / VtV) End Function '___________________________________________________________________________________________________________________________________________________________________________________________ '___________________________________________________________________________________________________________________________________________________________________________________________#End Region