| 일 | 월 | 화 | 수 | 목 | 금 | 토 |
|---|---|---|---|---|---|---|
| 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 |
- 프로그래밍
- CSS
- IT
- 위키
- 글이_점점_짧아지고_있습니다.
- 확률
- javascript
- 자바스크립트
- 나무위키
- C
- Java
- C++
- 프로그래밍 언어
- 수학
- 오늘의 문제
- 물리학
- html
- matlab
- ┌─┐
- HTML5.
- Nemeth Braille
- C언어
- html5
- 선형대수학
- LaTeX
- 점자
- python
- 리브레위키
- Inkscape
- css3
- Today
- Total
SturdyCobble's Study Note
MATLAB으로 (Reduced) Row Echelon Form 구하기 본문
NOTICE : 독학인 만큼, 잘못된 내용이 있을 수 있습니다. 내용 지적은 언제나 환영합니다.
더 이상 이 블로그에는 글을 올리지는 않을 예정입니다. 그렇지만 댓글을 달아주시면 최대한 답변드리고자 노력하겠습니다.
MATLAB에는 rref라는 함수로 위 Reduced Row Echelon Form을 구할 수 있습니다. 이번 글에서는 개략적으로 이 함수의 동작원리를 살펴보겠습니다. 자세한 사항은 https://en.wikipedia.org/wiki/Gaussian_elimination#Pseudocode 의 코드를 참조하시면 편합니다. (rref함수의 코드는 type rref로 구할 수 있습니다.)
function A = getRef(A)
% 1
[R,C] = size(A);
i = 1;
j = 1;
while i <= R && j <= C
% 2
[p, k] = max(abs(A(i:R,j)));
k = k+i-1;
% 3
if p == 0
j = j + 1;
else
% 4
A([i k],j:C) = A([k i],j:C);
% 5
A(i,j:C) = A(i,j:C)./A(i,j);
% 6
for k = [i+1:R]
A(k,j:C) = A(k,j:C) - A(k,j).*A(i,j:C);
end
% 7
i = i + 1;
j = j + 1;
end
end
위 코드는 단순히 행사다리꼴 행렬만을 구합니다. 물론 행사다리꼴은 한 개로 정해지지는 않지만 말이죠.
주석을 참고로 한 줄씩 해석해보겠습니다.
위의 행렬이 있다고 생각해봅시다. 그러면, 첫번째로 행렬의 크기를 구합니다.
그리고는 한 열벡터를 가져와 절댓값이 최대인 곳을 찾습니다. 이때 절댓값이 클수록 나눗셈을 할 때 부동소수점 오차가 줄어드는 효과가 있습니다. 또한, 절댓값의 최대를 찾음으로써 모든 값이 0인지 확인할 수 있습니다.
그리고 곧이어 k = k + i - 1로 인덱싱한 작은 행렬의 열벡터에서 구한 인덱스를 전체 행렬 기준으로 바꾸어줍니다.
만약 모든 값이 0이면, 다음 행에 대해서 조사해야 할 것입니다. 따라서 j값을 1 증가시킵니다.
0이 아닌 값이 있다면, 이를 맨 위로 올려줘야 할 것입니다. 이렇게 되면, 다음과 같이 행렬이 바뀝니다.
그다음 leading 1을 만들기 위해 나눠줘야 합니다. 그런데 왼쪽 위부터 연산을 진행하므로, i,j가 1이 아닌 상황을 생각해보면, 이미 다음과 같이 0이 존재할 것입니다. 따라서 그 오른쪽만 나눠줘도 충분합니다.
그 다음은 leading 1(i행 j열에 위치) 아래를 다 0으로 만들도록 적당히 수를 곱해 더해야 합니다.
그 다음으로는 leading 1을 새롭게 만들기 위해 i와 j를 1씩 증가시킵니다.
만약, reduced row echelon form을 구한다면, 아래와 같이 leading 1이 포함된 모든 행을 0으로 만들어야 할 것입니다.
for k = [1:i-1 i+1:R]
A(k,j:C) = A(k,j:C) - A(k,j).*A(i,j:C);
end(1:i-1 과 i+1:R의 두 행벡터를 결합하여 leading 1을 제외한 다른 열에 있는 요소를 0으로 만들었습니다.)
다음 코드로 reduce row echelon form을 찾으면서 그 과정을 살펴보겠습니다.
function A = getRref(A)
[R,C] = size(A);
i = 1;
j = 1;
while i <= R && j <= C
[p, k] = max(abs(A(i:R,j)));
k = k+i-1;
if p == 0
j = j + 1;
else
A([i k],j:C) = A([k i],j:C);
A
A(i,j:C) = A(i,j:C)./A(i,j);
for k = [1:i-1 i+1:R]
A(k,j:C) = A(k,j:C) - A(k,j).*A(i,j:C);
A
end
i = i + 1;
j = j + 1;
end
end
>> A = [ ...
1 2 3 4 5
7 8 9 10 11
11 12 13 14 15
15 16 17 18 19 ];
>>getRref(A)
A =
15 16 17 18 19
7 8 9 10 11
11 12 13 14 15
1 2 3 4 5
A =
1.0000 1.0667 1.1333 1.2000 1.2667
0 0.5333 1.0667 1.6000 2.1333
11.0000 12.0000 13.0000 14.0000 15.0000
1.0000 2.0000 3.0000 4.0000 5.0000
A =
1.0000 1.0667 1.1333 1.2000 1.2667
0 0.5333 1.0667 1.6000 2.1333
0 0.2667 0.5333 0.8000 1.0667
1.0000 2.0000 3.0000 4.0000 5.0000
A =
1.0000 1.0667 1.1333 1.2000 1.2667
0 0.5333 1.0667 1.6000 2.1333
0 0.2667 0.5333 0.8000 1.0667
0 0.9333 1.8667 2.8000 3.7333
A =
1.0000 1.0667 1.1333 1.2000 1.2667
0 0.9333 1.8667 2.8000 3.7333
0 0.2667 0.5333 0.8000 1.0667
0 0.5333 1.0667 1.6000 2.1333
A =
1.0000 0 -1.0000 -2.0000 -3.0000
0 1.0000 2.0000 3.0000 4.0000
0 0.2667 0.5333 0.8000 1.0667
0 0.5333 1.0667 1.6000 2.1333
A =
1.0000 0 -1.0000 -2.0000 -3.0000
0 1.0000 2.0000 3.0000 4.0000
0 0 -0.0000 -0.0000 -0.0000
0 0.5333 1.0667 1.6000 2.1333
A =
1.0000 0 -1.0000 -2.0000 -3.0000
0 1.0000 2.0000 3.0000 4.0000
0 0 -0.0000 -0.0000 -0.0000
0 0 0 0.0000 0
A =
1.0000 0 -1.0000 -2.0000 -3.0000
0 1.0000 2.0000 3.0000 4.0000
0 0 -0.0000 -0.0000 -0.0000
0 0 0 0.0000 0
A =
1.0000 0 0 -1.0625 -1.0000
0 1.0000 2.0000 3.0000 4.0000
0 0 1.0000 0.9375 2.0000
0 0 0 0.0000 0
A =
1.0000 0 0 -1.0625 -1.0000
0 1.0000 0 1.1250 0
0 0 1.0000 0.9375 2.0000
0 0 0 0.0000 0
A =
1.0000 0 0 -1.0625 -1.0000
0 1.0000 0 1.1250 0
0 0 1.0000 0.9375 2.0000
0 0 0 0.0000 0
A =
1.0000 0 0 -1.0625 -1.0000
0 1.0000 0 1.1250 0
0 0 1.0000 0.9375 2.0000
0 0 0 0.0000 0
A =
1.0000 0 0 0 -1.0000
0 1.0000 0 1.1250 0
0 0 1.0000 0.9375 2.0000
0 0 0 1.0000 0
A =
1.0000 0 0 0 -1.0000
0 1.0000 0 0 0
0 0 1.0000 0.9375 2.0000
0 0 0 1.0000 0
A =
1 0 0 0 -1
0 1 0 0 0
0 0 1 0 2
0 0 0 1 0
ans =
1 0 0 0 -1
0 1 0 0 0
0 0 1 0 2
0 0 0 1 0
위와 같이 깔끔하게 결과를 얻었습니다. 그러나 여기서 문제는 위 결과가 답이 아니라는 점입니다.
>> A = [ ...
1 2 3 4 5
7 8 9 10 11
11 12 13 14 15
15 16 17 18 19 ];
>> rref(A)
ans =
1 0 -1 -2 -3
0 1 2 3 4
0 0 0 0 0
0 0 0 0 0
여기서 문제는 부동소수점 숫자의 한계에서 옵니다. 정밀도에 한계가 있어 계산 결과가 0이 나와야 하는데, 0에 가까운 값이 구해져 0이 아니라 생각한 것이죠. 사실 내장 함수로도 정확히 0만 0으로 생각하게 하면 비슷한 결과가 나옵니다.
>> rref(A,0)
ans =
1 0 0 0 -1
0 1 0 0 0
0 0 1 0 2
0 0 0 1 0
이는 다음과 같이 정밀도를 지정함으로서 해결할 수 있습니다.
function A = ref01(A)
[R,C] = size(A);
i = 1;
j = 1;
while i <= R && j <= C
[p, k] = max(abs(A(i:R,j)));
k = k+i-1;
if p <= 0.0001
j = j + 1;
else
A([i k],j:C) = A([k i],j:C);
A
A(i,j:C) = A(i,j:C)./A(i,j);
for k = [1:i-1 i+1:R]
A(k,j:C) = A(k,j:C) - A(k,j).*A(i,j:C);
A
end
i = i + 1;
j = j + 1;
end
end
ans =
1.0000 0 -1.0000 -2.0000 -3.0000
0 1.0000 2.0000 3.0000 4.0000
0 0 -0.0000 -0.0000 -0.0000
0 0 0 0.0000 0
그래도 여러모로 부족한 부분이 보입니다. 이는 분모와 분자를 따로 생각한다던지 하는 방법으로 개선할 수 있습니다.
정밀도는 상황에 따라 다르지만, 내장 함수인 rref는 다음과 같이 정밀도를 정하고 있습니다.
tol = max(m,n)*eps(class(A))*norm(A,inf);
'프로그래밍 > MATLAB' 카테고리의 다른 글
| LDU Decomposition By MATLAB (0) | 2019.10.06 |
|---|---|
| [MATLAB] 10. 기초 연산 & 특수문자 (0) | 2019.09.08 |
| [MATLAB] 09. 데이터 처리와 파일 입출력 (0) | 2019.09.02 |
| [MATLAB] 08. 함수 (0) | 2019.09.01 |
| MATLAB으로 벡터장 표시하기 (0) | 2019.08.30 |
