盒子
盒子

矩阵求逆（sweep operator & Gauss Jordan elimination）

2015.10.08 nichunquan machine learning  热度 ℃

矩阵求逆我们一般可以使用Gauss—Jordan消去法，在此基础上，为了节省存储空间又可以使用sweep operator

####1. Gauss—Jordan elimination

对于矩阵A，将其扩展为[A|I]

P1

只进行每行的扩大缩小（乘除），行和行之间的加减

2

最后目的是将左矩阵变为单位矩阵，则右矩阵就是原矩阵的逆矩阵。这里我们会发现如图右矩阵每列一旦变成基向量，那么下面的变换这里就不再变化
这里的空间就被浪费了。可以使用sweep operator来解决。

####2. sweep operator

这里我们为了节省空间，可以将右矩阵的第一列存储到左矩阵的第一列，然后就可以不用存储右矩阵了。这样后面一次变换可以在前一次变换的基础上得出，

5

想要变回原矩阵，再次使用这个变化就行了。

sweep变换还有另外一种形式，只是求得的矩阵是逆矩阵的负值，这种变换叫做正sweep变换

3

把它变回来可以使用sweep变换的逆变换，

4

通过正sweep变换也可以求得矩阵的行列式的值，计算方法就是sweep变换的各步枢轴元素的乘积。

7

####3. python 和 R 代码

python:

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
import numpy as np
from numpy.linalg import matrix_rank
from numpy.linalg import det
def mySweep(A, m):
   	if A.shape[0] != A.shape[1]:
 	       return ('input matrix is not a squared matrix')
 	     elif A.shape[0] != matrix_rank(A):
   	    return ('input matrix is not invertible')
    	 else:
   	    B = np.copy(A)
   	    n = B.shape[0]
   	    for k in range(m):
   	        for i in range(n):
   	            for j in range(n):
   	                if i != k and j != k:
   	                    B[i,j] = B[i,j]- B[i,k]*B[k,j]/B[k,k]
   	        for i in range(n):
   	            if i != k:
   	                B[i,k] = B[i,k]/B[k,k]
   	        for j in range(n):
   	            if j != k:
   	                B[k,j] = B[k,j]/B[k,k]
   	        B[k,k] = - 1/B[k,k]
   	    return B
def mySweepRev(A, m):
   	if A.shape[0] != A.shape[1]:
   	    return ('input matrix is not a squared matrix')
   	elif A.shape[0] != matrix_rank(A):
   	    return ('input matrix is not invertible')
   	else:
   	    B = np.copy(A)
   	    n = B.shape[0]
   	    for k in range(m):
   	        for i in range(n):
   	            for j in range(n):
                   if i != k and j != k:
                       B[i,j] = B[i,j]- B[i,k]*B[k,j]/B[k,k]
           for i in range(n):
               if i != k:
                   B[i,k] = - B[i,k]/B[k,k]
           for j in range(n):
               if j != k:
                   B[k,j] = - B[k,j]/B[k,k]
           B[k,k] = - 1/B[k,k]
       return B
def myGaussJordan(A, m):
   	if A.shape[0] != A.shape[1]:
   	    return ('input matrix is not a squared matrix')
   	elif A.shape[0] != matrix_rank(A):
   	    return ('input matrix is not invertible')
   	else:
   	    n = A.shape[0]
   	    B = np.hstack((A, np.identity(n)))
   	    for k in range(m):
   	        a = B[k, k]
   	        for j in range(n * 2):
   	            B[k, j] = B[k, j] / a
   	        for i in range(n):
   	            if i != k:
   	                a = B[i, k]
   	                for j in range(n * 2):
   	                    B[i, j] = B[i, j] - B[k, j] * a;
   	    return B[:, n:2 * n]
def myGaussJordanVec(A, m):
   	if A.shape[0] != A.shape[1]:
   	    return ('input matrix is not a squared matrix')
   	elif A.shape[0] != matrix_rank(A):
   	    return ('input matrix is not invertible')
   	else:
   	    n = A.shape[0]
   	    B = np.hstack((A, np.identity(n)))
   	    for k in range(m):
   	        B[k, :] = B[k,] / B[k, k]
   	        for i in range(n):
   	            if i != k:
   	                B[i,] = B[i,] - B[k,] * B[i, k]
   	    return B[:, n:2 * n]
def mySweepDet(A, m):
   	if A.shape[0] != A.shape[1]:
   	    return ('input matrix is not a squared matrix')
   	elif A.shape[0] != matrix_rank(A):
   	    return ('input matrix is not invertible')
   	else:
   	    B = np.copy(A)
   	    D = 1
   	    n = B.shape[0]
   	    for k in range(m):
   	        for i in range(n):
   	            for j in range(n):
   	                if i != k and j != k:
   	                    B[i,j] = B[i,j]- B[i,k]*B[k,j]/B[k,k]
   	        for i in range(n):
   	            if i != k:
   	                B[i,k] = - B[i,k]/B[k,k]
   	        for j in range(n):
   	            if j != k:
   	                B[k,j] =  B[k,j]/B[k,k]
   	        D = D * B[k,k]
   	        B[k,k] =  1/B[k,k]
   	    return D
A = np.array([[1,2,3],[7,11,13],[17,21,23]], dtype=float).T
print('sweep operator:\n',mySweep(A, 3))
print('reverse sweep operator:\n',mySweepRev(mySweep(A, 3),3))
print('GaussJordan:\n',myGaussJordan(A, 3))
print('GaussJordan vectorized version:\n',myGaussJordanVec(A, 3))
print('computing the determinant of a matrix using the det function:',det(A))
print('computing the determinant of a matrix using the sweep operator:',mySweepDet(A, 3))

R code:

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
mySweep <- function(A, m)
{	
  if(nrow(A)==ncol(A)&&nrow(A)==qr(A)$rank)
	 {
   	n <- dim(A)[1]
   	
   	for (k in 1:m) 
   	{
   	  for (i in 1:n)     
   	    for (j in 1:n)   
   	      if (i!=k  & j!=k)     
   	        A[i,j] <- A[i,j] - A[i,k]*A[k,j]/A[k,k]    
   	  
   	  for (i in 1:n) 
   	    if (i!=k) 
   	      A[i,k] <- A[i,k]/A[k,k]  
   	  
   	  for (j in 1:n) 
   	    if (j!=k) 
   	      A[k,j] <- A[k,j]/A[k,k]  
   	  
   	  A[k,k] <- - 1/A[k,k] 
   	}
   	return(A)
 	 }	
	 else
    return('input matrix is not a squared matrix or not invertible' )
}
mySweepRev <- function(A, m)
{
	 if(nrow(A)==ncol(A)&&nrow(A)==qr(A)$rank)
	 {
   	n <- dim(A)[1]
   	
   	for (k in 1:m) 
   	{
   	  for (i in 1:n)     
   	    for (j in 1:n)   
   	      if (i!=k  & j!=k)     
   	        A[i,j] <- A[i,j] - A[i,k]*A[k,j]/A[k,k]    
   	  
   	  for (i in 1:n) 
   	    if (i!=k) 
   	      A[i,k] <- - A[i,k]/A[k,k]  
   	  
   	  for (j in 1:n) 
   	    if (j!=k) 
   	      A[k,j] <- - A[k,j]/A[k,k]  
   	  
   	  A[k,k] <- - 1/A[k,k] 
     }
   	return(A)
    }	
 	 else
    return('input matrix is not a squared matrix or not invertible' )
}
myGaussJordan <- function(A, m)
{
	 if(nrow(A)==ncol(A)&&nrow(A)==qr(A)$rank)
	 {
   	n <- dim(A)[1]
   	B <- cbind(A, diag(rep(1, n)))
   	
   	for (k in 1:m) 
   	{
   	  a <- B[k, k]
   	  for (j in 1:(n*2))     
   	    B[k, j] <- B[k, j]/a
   	  for (i in 1:n)
   	    if (i != k)
   	    {
         a <- B[i, k]
   	      for (j in 1:(n*2))
   	        B[i, j] <- B[i, j] - B[k, j]*a; 
   	    }    
   	}
   	return(B)
	 }
	 else
   	return('input matrix is not a squared matrix or not invertible' )
 
}
myGaussJordanVec <- function(A, m)
{
 		if(nrow(A)==ncol(A)&&nrow(A)==qr(A)$rank)
   {
   	n <- dim(A)[1]
   	B <- cbind(A, diag(rep(1, n)))
   	
   	for (k in 1:m) 
   	{
   	  B[k, ] <- B[k, ]/B[k, k]
   	  for (i in 1:n)
       if (i != k)
   	      B[i, ] <- B[i, ] - B[k, ]*B[i, k];   
   	}
   	return(B)
   }	
	 else
   	return('input matrix is not a squared matrix or not invertible' )
   }
mySweepDet <- function(A, m)
{
  if(nrow(A)==ncol(A)&&nrow(A)==qr(A)$rank)
  {
   	
   	n <- dim(A)[1]
   	D <- 1
   	for (k in 1:m) 
   	{
   	  for (i in 1:n)     
   	    for (j in 1:n)   
   	      if (i!=k  & j!=k)     
   	        A[i,j] <- A[i,j] - A[i,k]*A[k,j]/A[k,k]    
     
   	  for (i in 1:n) 
   	    if (i!=k) 
   	      A[i,k] <- A[i,k]/A[k,k]  
   	  
   	  for (j in 1:n) 
   	    if (j!=k) 
   	      A[k,j] <- A[k,j]/A[k,k]  
   	  D <- D * A[k,k]
   	  A[k,k] <- - 1/A[k,k] 
   	}
   	return(D)
  }	
  else
   	return('input matrix is not a squared matrix or not invertible' )
}
A = matrix(c(1,2,3,7,11,13,17,21,23), 3,3)
solve(A)
mySweep(A,3)
mySweepRev(mySweep(A,3),3)
myGaussJordan(A,3)
myGaussJordanVec(A,3)
mySweepDet(A,3)