C++

Forum

 
matrix inverse

moctar (9)
i am trying to write a program to get the inverse of any type of square matrix using strassen based matrix inversion.
but what i don't understand is how to represent the sub-block matrix in the form A^-1, B^-1 etc...

since i have to apply it to different matrixes i can't really use something like
A^-1 = 1/determinant. so how to implement the blockkwise inverse.
a code would really be helpful.
lastchance (6980)
What do you want to do with this matrix inverse?

Note that it's actually surprisingly rare to require a matrix inverse. It's usually easier to QR or LU factorise, because those factorised systems can be inverted quickly as and when necessary.
moctar (9)
i mean not a program but a function, to get the inverse of a matrix using strassen matrix inversion method given a matrix as input. It is explained on this link but i don't know how to implement it.

http://tesla.pmf.ni.ac.rs/people/DeXteR/old/Papers/RapidLU.pdf
lastchance (6980)
You haven't answered my question.

What are you intending to do with the matrix inverse once you've got it?
moctar (9)
i want to multiply the matrix inverse with the matrix using strassen multiplication to get the identity matrix
lastchance (6980)
i want to multiply the matrix inverse with the matrix using strassen multiplication to get the identity matrix

Well, I can think of easier ways of getting the identity matrix.


Anyway, I tried coding from the link that you provided (see the Algorithm 2.1 on page 4, and the equations it refers to.)

It appears to work ...

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
 
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <cmath>
using namespace std;

const double SMALL = 1.0E-30;          // used to stop divide-by-zero
const double NEARZERO = 1.0e-10;       // helps in printing

using vec    = vector<double>;         // vector
using matrix = vector<vec>;            // matrix

// Function prototypes
void print( const string &title, const matrix &A );
matrix matmul( const matrix &A, const matrix &B );
matrix subtract( const matrix &A, const matrix &B );
matrix oppsign( matrix A );
matrix subMatrix( const matrix &A, int i1, int i2, int j1, int j2 );
matrix assembly( const matrix &A11, const matrix &A12, const matrix &A21, const matrix &A22 );
matrix inverse( const matrix &A );

//======================================================================

void print( const string &title, const matrix &A )
{
   if ( title != "" ) cout << title << '\n';

   for ( auto &row : A )
   {
      for ( auto x : row ) cout << setw( 15 ) << ( abs( x ) < NEARZERO ? 0.0 : x );
      cout << '\n';
   }
}

//======================================================================

matrix matmul( const matrix &A, const matrix &B )          // Matrix times matrix
{
   int rowsA = A.size(),   colsA = A[0].size();
   int rowsB = B.size(),   colsB = B[0].size();
   assert( colsA == rowsB );

   matrix C( rowsA, vec( colsB, 0.0 ) );
   for ( int i = 0; i < rowsA; i++ )
   {
      for ( int j = 0; j < colsB; j++ )
      {
         for ( int k = 0; k < colsA; k++ ) C[i][j] += A[i][k] * B[k][j];
      }
   }
   return C;
}

//======================================================================

matrix subtract( const matrix &A, const matrix &B )        // Subtract matrices
{
   int rows = A.size(),   cols = A[0].size();
   assert( rows == B.size() && cols == B[0].size() );

   matrix result( rows, vec( cols ) );
   for ( int i = 0; i < rows; i++ )
   {
      for ( int j = 0; j < cols; j++ ) result[i][j] = A[i][j] - B[i][j];
   }
   return result;
}

//======================================================================

matrix oppsign( matrix A )                                  // Minus matrix
{
   for ( auto &row : A )
   {
      for ( auto &e : row ) e = -e;
   }
   return A;
}

//======================================================================

matrix subMatrix( const matrix &A, int i1, int i2, int j1, int j2 )
{
   int rows = i2 - i1 + 1, cols = j2 - j1 + 1;
   matrix result( rows, vec( cols ) );
   for ( int i = i1, r = 0; i <= i2; i++, r++ )
   {
      auto it1 = A[i].begin() + j1, it2 = A[i].begin() + j2 + 1;
      copy( it1, it2, result[r].begin() );
   }
   return result;
}

//======================================================================

matrix assembly( const matrix &A11, const matrix &A12, const matrix &A21, const matrix &A22 )
{
   int k = A11.size();           
   int n = k + A22.size();
   matrix result( n, vec( n ) );

   for ( int i = 0; i < k; i++ )
   {
      copy( A11[i].begin(), A11[i].end(), result[i].begin()     );
      copy( A12[i].begin(), A12[i].end(), result[i].begin() + k );
   }

   for ( int i = k; i < n; i++ )
   {
      copy( A21[i-k].begin(), A21[i-k].end(), result[i].begin()     );
      copy( A22[i-k].begin(), A22[i-k].end(), result[i].begin() + k );
   }

   return result;
}

//======================================================================

matrix inverse( const matrix &A )
{
   int n = A.size();
   if ( n == 1 ) 
   {
      double value = A[0][0];
      if ( abs( value ) < SMALL )
      {
         cerr << "Non-invertible. Giving up.\n";
         exit( 0 );
      }
      return matrix( 1, vec( 1, 1.0 / value ) );
   }

   // Partition into four
   int k = n / 2;
   matrix A11 = subMatrix( A, 0, k - 1, 0, k - 1 );
   matrix A12 = subMatrix( A, 0, k - 1, k, n - 1 );
   matrix A21 = subMatrix( A, k, n - 1, 0, k - 1 );
   matrix A22 = subMatrix( A, k, n - 1, k, n - 1 );

   // Strassen steps
   matrix R1  = inverse( A11 );
   matrix R2  = matmul( A21, R1 );
   matrix R3  = matmul( R1, A12 );
   matrix R4  = matmul( A21, R3 );
   matrix R5  = subtract( R4, A22 );
   matrix R6  = inverse( R5 );
   matrix X12 = matmul( R3, R6 );
   matrix X21 = matmul( R6, R2 );
   matrix R7  = matmul( R3, X21 );
   matrix X11 = subtract( R1, R7 );
   matrix X22 = oppsign( R6 );

   return assembly( X11, X12, X21, X22);
}

//====================================================================== 

int main()
{
   // Data
   matrix A = { { 1, 2, 3, 4, 5 }, { 9, 88, 7, 63, 5 }, { 1, 3, 15, 32, 1 }, { 2, 4, 22, 4, 222 }, { 3, 5, 9, 0, 0 } };
   print( "A:", A );
   matrix B = inverse( A );
   print( "\nB:", B );
   print( "\nCheck AB=I:", matmul( A, B ) );
}


A:
              1              2              3              4              5
              9             88              7             63              5
              1              3             15             32              1
              2              4             22              4            222
              3              5              9              0              0

B:
        2.34541     -0.0288739      -0.229938     -0.0511385       -0.25111
       -0.36173      0.0150113       0.014695     0.00774279      0.0654828
      -0.580842     0.00128504      0.0684823      0.0127446       0.158435
       0.231676    -0.00110368     0.00511536    -0.00521612     -0.0721421
      0.0387744   -0.000117808    -0.00507196     0.00365671     -0.0133185

Check AB=I:
              1              0              0              0              0
              0              1              0              0              0
              0              0              1              0              0
              0              0              0              1              0
              0              0              0              0              1


Last edited on
moctar (9)
Thank you! i will have a look into it.
Topic archived. No new replies allowed.