70강에서 소거의 결과를 저장해 재사용했습니다. 상수항이 바뀌어도 분해는 그대로이므로, 한 번 하고 여러 번 씁니다.
80강의 그람슈미트도 같은 성격입니다. 직교화를 한 번 해 두면 좌표 계산과 정사영이 모두 싸집니다. 그러면 그 결과를 어떻게 저장해야 합니까.
답은 80강 심화 1에 있었습니다. 그람슈미트를 실행하면서 나온 계수를 모으면 위삼각행렬이 되고
이 성립합니다. 이 강의에서 그 분해를 정식으로 다룹니다.
그리고 왜 이 분해가 최소제곱의 표준 방법인지를 밝힙니다. 82강에서 배울 정규방정식은 조건수를 제곱하는데, QR은 그러지 않습니다. 검산에서 그 차이가 **오차 과 **으로 나타납니다.
문제. 의 열이 , , 입니다.
(1) 80강의 그람슈미트를 실행하며 나온 계수를 모두 적으세요.
(2) 그 계수들로 행렬 을 만드세요.
(3) 을 계산하세요.
생각의 실마리. 그람슈미트에서 두 종류의 수가 나옵니다. 빼는 데 쓴 내적 값과 정규화할 때 나눈 길이입니다. 둘 다 저장합니다.
풀이. (1) 80강 문제 2의 계산에서
(2) 위치를 정합니다. ** 성분을 **로 두면
검산에서 확인되며 위삼각이고 대각이 양수입니다.
(3) 입니다. 검산에서 확인됩니다.
이 문제에서 배우는 것: QR 분해.
QR 분해. 열이 일차독립인 (, )에 대해
인 (, 정규직교 열)와 위삼각 (, 대각 양수)가 유일하게 존재합니다.
의 성분이 무엇인지 봅니다.
이므로 입니다. 검산에서 가 과 같습니다.
위삼각인 이유는 80강 문제 2의 부분공간 보존입니다. 가 로만 만들어지므로, 이면 입니다.
대각 성분의 뜻도 있습니다.
정규화하기 전 벡터의 길이이며, 가 앞의 것들과 얼마나 독립적인지를 잽니다. 작으면 거의 종속입니다.
행렬식과의 관계도 나옵니다. 가 직교이므로 이고
검산에서 둘 다 입니다. 79강의 부피 배율이 여기서 다시 계산됩니다.
70강의 LU와 비교합니다.
| 분해 | 만드는 법 | 또는 | 비용 |
|---|---|---|---|
| 소거 | 아래삼각 | ||
| 직교화 | 직교 |
QR이 두 배 비싼데 훨씬 안정합니다. 문제 5에서 그 값어치를 봅니다.
바로 확인 1.
확인 1-1. 의 성분을 쓰세요.
답. 이며 입니다.
확인 1-2. 이 위삼각인 이유를 쓰세요.
답. 가 처음 개의 로만 만들어지기 때문입니다.
확인 1-3. 가 무엇을 뜻합니까?
답. 정규화 전 벡터의 길이이며 앞의 것들과의 독립성 정도입니다.
문제. 가 이고 열이 과 입니다.
(1) 그람슈미트로 얻는 와 의 크기를 쓰세요.
(2) 를 정사각으로 만들 수 있습니까?
(3) 두 형태를 비교하세요.
생각의 실마리. 그람슈미트는 열마다 하나씩 를 만들므로 열 개수만큼만 나옵니다. 그런데 의 완전한 기저는 네 개가 필요합니다.
풀이. (1) 검산에서 가 이고 이 입니다.
(2) 됩니다. 의 정규직교기저 두 개를 더 붙이면 가 됩니다.
검산에서 mode='complete'로 가 이고 입니다. 이때 은 이며 아래 두 행이 입니다.
(3) 정리합니다.
이 문제에서 배우는 것: 두 가지 QR.
축소형과 완전형. 가 ()일 때
형태 관계 축소형 Q^{\top}Q=I_ 완전형 Q^{\top}Q=QQ^{\top}=I_
축소형의 는 정사각이 아니므로 직교행렬이 아닙니다. 이지만 입니다.
실제로 은 80강 문제 3의 정사영 행렬입니다.
74강의 네 부분공간이 완전형에서 한눈에 드러납니다.
| 의 열 | 생성하는 공간 |
|---|---|
| 앞 개 | |
| 나머지 개 |
완전형 QR이 열공간과 좌영공간의 정규직교기저를 동시에 줍니다. 74강 심화 6에서 "실무에서는 분해로 기저를 구한다"고 한 것이 이것입니다.
의 아래쪽이 인 것도 자연스럽습니다. 의 열이 앞 개 로만 만들어지므로, 뒤의 와의 내적이 입니다.
실무에서는 대개 축소형을 씁니다. 필요한 것이 열공간이고, 이 크면 행렬을 저장할 이유가 없습니다.
바로 확인 2.
확인 2-1. 축소형 QR에서 의 크기를 쓰세요.
답. 와 같은 입니다.
확인 2-2. 축소형 가 직교행렬입니까?
답. 아닙니다. 이지만 입니다.
확인 2-3. 이 무엇입니까?
답. 열공간으로의 정사영 행렬입니다.
문제. 문제 2의 와 를 봅니다.
(1) 에 해가 있습니까?
(2) QR로 최소제곱 해를 구하세요.
(3) 잔차를 확인하세요.
생각의 실마리. 가 이므로 방정식이 넷이고 미지수가 둘입니다. 68강 심화 5의 과결정계이며 대개 해가 없습니다.
풀이. (1) 없습니다. 이므로 열공간이 의 이차원 부분공간이고, 가 거기 있을 이유가 없습니다.
(2) 를 최소로 하는 를 찾습니다. 을 넣으면
인데, 80강에서 직교변환이 노름을 보존하므로 을 곱해도 됩니다. 자세한 유도는 문제 4에서 하고, 결과는
입니다. 이 위삼각이므로 후진대입 한 번이면 끝입니다.
검산에서 이고, 정규방정식으로 푼 값과 같습니다.
(3) 잔차 에 대해 입니다. 검산에서 정확히 영벡터이고 입니다.
이 문제에서 배우는 것: QR로 푸는 최소제곱.
절차.
- (축소형)
- \mathbf{c}=Q^{\top}\mathbf
- 를 후진대입으로 풉니다
분해가 비싸고 나머지는 쌉니다. 70강의 구조와 같습니다.
| 단계 | 비용 |
|---|---|
| 분해 | 2mn^ |
| Q^{\top}\mathbf | |
| 후진대입 | n^ |
가 여러 개면 분해를 재사용합니다. 같은 설계행렬에 여러 반응변수를 회귀할 때 쓰는 방법입니다.
잔차가 인 것은 잔차가 좌영공간에 있다는 뜻이며, 74강 심화 4에서 예고한 내용입니다.
82강에서 이 조건이 정규방정식이 됩니다. 여기서는 QR로 푼 결과가 그 조건을 만족함을 확인했습니다.
바로 확인 3.
확인 3-1. QR로 최소제곱을 푸는 식을 쓰세요.
답. 입니다.
확인 3-2. 이 위삼각인 것이 왜 유리합니까?
답. 후진대입 한 번으로 풀리기 때문입니다.
확인 3-3. 잔차가 만족하는 조건을 쓰세요.
답. 이며 좌영공간에 있습니다.
문제. 를 QR로 유도합니다.
(1) 완전형 을 쓰고 을 곱하세요.
(2) 노름이 어떻게 갈라집니까?
(3) 최소가 되는 조건을 쓰세요.
생각의 실마리. 직교변환이 노름을 보존하므로 을 곱해도 최소화 문제가 바뀌지 않습니다. 그러면 문제가 훨씬 단순해집니다.
풀이. (1) 완전형에서 가 직교행렬이고 이 입니다. 75강 문제 3에서
(2) 의 아래 개 행이 이므로 위아래로 가릅니다.
그러면
(3) 둘째 항은 와 무관합니다. 손댈 수 없습니다. 첫째 항을 으로 만들면 최소이며
입니다. 이 가역(대각이 양수)이므로 해가 유일합니다.
이 문제에서 배우는 것: 직교변환이 문제를 단순화합니다.
핵심. 직교변환은 노름을 보존하므로 최소화 문제를 바꾸지 않고 형태만 단순하게 만듭니다.
이 논법이 이 과목에서 반복됩니다.
| 문제 | 직교변환으로 |
|---|---|
| 최소제곱 (81강) | 삼각계로 |
| 대칭행렬 (86강) | 대각으로 |
| 일반 행렬 (88강) | 대각으로 (기저 둘) |
**"보존하면서 단순화"**가 직교변환의 값어치입니다. 일반 가역변환은 단순화는 하지만 노름을 바꾸므로, 최소화 문제에 쓸 수 없습니다.
기하적으로도 읽힙니다. 74강 문제 4의 그림에서
이고 는 언제나 안에 있으므로 을 만들 수 없습니다. ****이며, 그것이 피할 수 없는 최소 오차입니다.
가 의 뒤쪽 열, 즉 좌영공간 기저와의 내적이므로 정확히 맞습니다.
바로 확인 4.
확인 4-1. 을 곱해도 되는 이유를 쓰세요.
답. 직교변환이 노름을 보존하기 때문입니다.
확인 4-2. 최소 오차가 무엇입니까?
답. 이며 의 좌영공간 성분의 크기입니다.
확인 4-3. 왜 를 줄일 수 없습니까?
답. 와 무관한 항이기 때문입니다.
문제. 이 작을 때 을 봅니다. 참해가 이 되도록 로 둡니다.
(1) 정규방정식으로 풀고 오차를 재세요.
(2) QR로 풀고 오차를 재세요.
(3) 비교하세요.
생각의 실마리. 두 열이 거의 같습니다. 를 만들면 그 유사성이 제곱으로 증폭됩니다.
풀이. 검산에서 세 경우를 비교합니다.
| 정규방정식 오차 | QR 오차 | ||
|---|---|---|---|
| 10^ | 1.414\times10^ | 1.014\times10^ | 2.220\times10^ |
| 10^ | 1.414\times10^ | 1.171\times10^ | |
| 10^ | 1.414\times10^ | 1.131\times10^ | 2.220\times10^ |
정규방정식의 오차가 에 비례해 커집니다. 이 배 작아질 때마다 오차가 약 배가 됩니다.
QR은 기계 정밀도에 머뭅니다.
이 문제에서 배우는 것: 조건수를 제곱하지 마세요.
핵심. 정규방정식은 이므로 조건수를 제곱합니다. QR은 를 유지합니다.
왜 제곱이 되는가. 88강에서 특이값으로 설명되지만 미리 말하면, 의 특이값이 일 때 의 고유값이 입니다. 그러면
입니다.
실용적 결론입니다.
배정밀도에서 이면 정규방정식은 유효숫자가 남지 않고 QR은 여덟 자리가 남습니다.
방법을 정리합니다.
| 방법 | 비용 | 안정성 | 쓰는 곳 |
|---|---|---|---|
| 정규방정식 + 콜레스키 | \kappa^ | 조건이 좋고 빠름이 중요할 때 | |
| QR | 2mn^ | 표준 | |
| SVD | 계수 부족까지 다룰 때 |
검산에서 세 크기의 비용을 비교합니다.
| QR | 정규방정식 | ||
|---|---|---|---|
| 1.933\times10^ | 1.033\times10^ | ||
| 4.167\times10^ | 2.917\times10^ | ||
| 3.947\times10^ | 2.027\times10^ |
QR이 약 두 배 비쌉니다. 그 대가로 조건수를 제곱하지 않습니다.
대부분의 통계 소프트웨어가 QR을 씁니다. R의 lm, 파이썬의 lstsq가 그렇습니다. 82강에서 정규방정식을 배우지만, 그것은 이해를 위한 것이고 계산은 QR로 합니다.
69강 문제 5, 70강 문제 3, 78강 심화 6에 이어 다섯 번째로 같은 교훈입니다.
바로 확인 5.
확인 5-1. 정규방정식의 조건수를 쓰세요.
답. 입니다.
확인 5-2. QR의 조건수를 쓰세요.
답. 입니다.
확인 5-3. 두 방법의 비용 차이를 쓰세요.
답. QR이 약 두 배 비쌉니다.
| 개념 | 내용 |
|---|---|
| QR 분해 | , 정규직교 열, 위삼각 |
| 의 성분 | , 즉 |
| 대각 성분 | |
| 행렬식 | \lvert\det A\rvert=\prod r_ |
| 축소형 | 가 , Q^{\top}Q=I_ |
| 완전형 | 가 직교행렬 |
| (축소형) | 열공간 정사영 |
| 최소제곱 | 절차 |
|---|---|
| 1 | |
| 2 | \mathbf{c}=Q^{\top}\mathbf |
| 3 | 후진대입 |
| 최소 오차 |
| 안정성 | 조건수 | 비용 |
|---|---|---|
| 정규방정식 | \kappa^ | |
| QR | 2mn^ |
| 자주 하는 실수 | 바로잡기 |
|---|---|
| 축소형 를 직교행렬이라 합니다 | 입니다 |
| 을 아래삼각으로 씁니다 | 위삼각입니다 |
| 정규방정식을 코드로 씁니다 | QR을 씁니다 |
| 완전형을 기본으로 씁니다 | 대개 축소형이면 충분합니다 |
문제 6. QR 분해의 정의를 쓰세요.
답. 이며 의 열이 정규직교이고 이 대각 양수인 위삼각입니다.
문제 7. 을 와 로 나타내세요.
답. 입니다.
문제 8. 이 위삼각인 이유를 쓰세요.
답. 가 처음 개의 로만 만들어지기 때문입니다.
문제 9. 가 무엇을 뜻합니까?
답. 정규화 전 벡터의 길이이며 앞의 것들과의 독립성 정도입니다.
문제 10. 를 로 쓰세요.
답. 입니다.
문제 11. 축소형과 완전형의 크기를 쓰세요.
답. 각각 과 입니다.
문제 12. 축소형에서 이 무엇입니까?
답. 열공간으로의 정사영 행렬입니다.
문제 13. 완전형 의 뒤쪽 열이 무엇을 생성합니까?
답. 좌영공간 입니다.
문제 14. QR로 최소제곱을 푸는 식을 쓰세요.
답. 입니다.
문제 15. 을 곱해도 되는 이유를 쓰세요.
답. 직교변환이 노름을 보존하기 때문입니다.
문제 16. 최소 오차가 무엇입니까?
답. 의 좌영공간 성분의 크기입니다.
문제 17. 정규방정식과 QR의 조건수를 비교하세요.
답. 과 입니다.
문제 18. 가 여러 개일 때 무엇을 재사용합니까?
답. 분해 을 재사용합니다.
심화 1. QR 분해의 유일성을 증명하세요.
풀이. 대각이 양수라는 조건을 붙이면 유일합니다.
라 하고 두 이 대각 양수인 위삼각이라 합니다.
왼쪽을 봅니다. 가 같은 열공간의 정규직교기저이므로, 이 이고
입니다. 가운데에서 가 열공간 정사영이고 의 열이 그 안에 있으므로 입니다. 따라서 이 직교행렬입니다.
오른쪽을 봅니다. 위삼각행렬의 역과 곱이 위삼각이므로 이 위삼각입니다.
직교이면서 위삼각이면 대각행렬입니다. 의 첫 열을 보면 이고, 첫 행이 다른 열과 직교하므로 나머지가 정해집니다. 귀납적으로 진행하면 대각뿐입니다.
대각 성분이 인데, 에서 대각이 이고 둘 다 양수여야 하므로 입니다.
**"직교이면서 삼각이면 대각"**이라는 보조 사실이 핵심이었습니다. 이것은 뒤에서도 쓰입니다.
70강의 LU와 비교하면 조건이 비슷합니다. LU도 의 대각을 로 고정해야 유일했습니다. 정규화 조건이 유일성을 만듭니다.
심화 2. QR 알고리즘으로 고유값을 구하는 원리를 밝히세요.
풀이. QR 분해의 가장 놀라운 응용입니다.
QR 알고리즘. 에서 시작해
를 반복하면 가 (대개) 삼각행렬에 수렴하고, 대각에 고유값이 나타납니다.
왜 고유값이 보존되는가. 이므로
직교 닮음입니다. 77강에서 닮음이 고유값을 보존한다고 했으므로, 모든 가 같은 고유값을 가집니다.
왜 수렴하는가는 더 깊은 이야기입니다. 거듭제곱법과 관련되며, 큰 고유값 방향이 점점 앞으로 밀려 나옵니다.
검산에서 확인합니다. 대칭행렬
에 번 반복하면 대각이
이고, eigvalsh의 결과와 소수 여덟 자리까지 일치합니다.
이 알고리즘이 실무의 표준입니다. 84강에서 고유값을 특성다항식의 근으로 정의하지만, 실제로는 다항식을 풀지 않습니다. 이유가 있습니다.
| 방법 | 문제 |
|---|---|
| 특성다항식의 근 | 계수 계산과 근 구하기가 모두 불안정 |
| QR 알고리즘 | 직교 변환만 쓰므로 안정 |
다항식의 근은 계수에 극도로 민감합니다. 윌킨슨의 예에서 차수 인 다항식의 계수를 만 바꿔도 근이 크게 움직입니다. 그래서 거꾸로 갑니다. 다항식의 근이 필요하면 그것을 특성다항식으로 갖는 동반행렬을 만들어 QR 알고리즘을 돌립니다.
실무 구현에는 개선이 여럿 들어갑니다.
| 기법 | 목적 |
|---|---|
| 헤센베르크 축소 | 반복 비용을 에서 으로 |
| 이동 | 수렴 속도를 높입니다 |
| 수축 | 수렴한 고유값을 떼어냅니다 |
이 개선까지 넣으면 고유값 계산이 약 입니다.
심화 3. 열 피벗팅 QR을 소개하고 계수 판정에 쓰는 법을 논하세요.
풀이. 문제 1에서 가 독립성의 정도라 했습니다. 그것을 이용합니다.
열 피벗팅 QR. 매 단계에서 남은 부분의 노름이 가장 큰 열을 다음으로 고릅니다.
여기서 는 열 순열행렬입니다.
그러면 의 대각이 감소하는 순서가 됩니다.
계수 판정에 쓸 수 있습니다. 가 문턱보다 작아지는 지점에서 자르면 그것이 수치적 계수입니다.
67강 문제 5의 부분 피벗팅과 같은 발상입니다. 큰 것을 먼저 처리해 작은 것으로 나누는 일을 피합니다.
용도를 정리합니다.
| 용도 | 내용 |
|---|---|
| 계수 판정 | 대각이 언제 작아지는가 |
| 열 선택 | 앞쪽 열들이 중요한 열 |
| 기저 추출 | 열공간 기저를 원래 열에서 |
| 최소제곱 (계수 부족) | 축퇴를 처리 |
둘째와 셋째 줄이 실무에서 유용합니다. 72강 심화 1에서 열공간 기저를 원래 열에서 골라야 한다고 했는데, 소거보다 이 방법이 안정합니다.
SVD와 비교하면 이렇습니다.
| 열 피벗팅 QR | SVD | |
|---|---|---|
| 계수 판정 | 대개 정확 | 가장 정확 |
| 비용 | 이지만 상수가 큽니다 | |
| 원래 열 유지 | 예 | 아닙니다 |
셋째 줄이 결정적인 경우가 있습니다. 특징 선택에서는 "어느 변수가 중요한가"를 원하는데, SVD의 특이벡터는 원래 변수의 조합이라 해석이 어렵습니다. 열 피벗팅 QR은 원래 열을 고르므로 답이 바로 해석됩니다.
심화 4. QR 분해를 갱신하는 방법을 논하세요.
풀이. 데이터가 하나씩 들어오는 상황을 봅니다. 에 행이 추가될 때마다 처음부터 분해하면 낭비입니다.
행 추가. 에 행 이 붙으면
오른쪽이 위삼각이 아니므로 기븐스 회전으로 아래 행을 없앱니다. 번의 회전이면 되므로 입니다.
이 크면 결정적인 차이입니다.
열 추가와 삭제도 비슷하게 처리됩니다.
| 연산 | 비용 |
|---|---|
| 행 추가 | |
| 행 삭제 | |
| 열 추가 | |
| 열 삭제 |
응용을 봅니다.
| 상황 | 갱신되는 것 |
|---|---|
| 온라인 회귀 | 데이터가 하나씩 들어옵니다 |
| 이동창 회귀 | 오래된 데이터를 버립니다 |
| 단계적 변수 선택 | 변수를 넣고 뺍니다 |
| 교차검증 | 표본 하나를 뺍니다 |
넷째 줄이 유용합니다. 하나 남기기 교차검증에서 번 재계산하는 대신 갱신하면 훨씬 쌉니다. 69강 심화 4의 셔먼-모리슨과 같은 발상이며, 분해를 갱신하는 편이 수치적으로 더 안전합니다.
주의도 있습니다. 갱신을 반복하면 오차가 쌓입니다. 일정 횟수마다 다시 분해하는 것이 보통이며, 70강 심화 6에서 한 이야기와 같습니다.
심화 5. 최소제곱의 잔차와 통계적 해석을 잇습니다.
풀이. 문제 3에서 잔차가 좌영공간에 있다고 했습니다. 통계에서 무슨 뜻인지 봅니다.
회귀 모형 에서 최소제곱 추정은
입니다. 80강의 정사영 분해 그대로입니다.
자유도가 여기서 나옵니다. 잔차가 사는 공간이 이고 그 차원이 이므로
**73강 심화 3에서 "통계의 자유도가 퇴화차수"**라고 한 것이 이것입니다. 분산 추정에서 이 아니라 로 나누는 이유입니다.
피타고라스 정리가 분산분석이 됩니다.
평균을 뺀 형태로 쓰면 총제곱합이 회귀제곱합과 잔차제곱합으로 갈라지며, 이 그 비율입니다.
64강 심화 3에서 평균 빼기가 방향 정사영이라 했는데, 여기서 그 조작이 다시 나옵니다.
QR이 통계 계산에 주는 이점도 있습니다.
| 필요한 것 | QR로 |
|---|---|
| \hat | R^{-1}Q^{\top}\mathbf |
| 잔차 | \mathbf{y}-Q Q^{\top}\mathbf |
| (X^{\top}X)^ | R^{-1}R^ |
| 레버리지 | 의 대각 |
셋째 줄이 표준오차 계산에 필요한데, 를 만들지 않고 만으로 얻습니다. 조건수를 제곱하지 않는 이점이 여기까지 이어집니다.
넷째 줄의 레버리지는 각 관측이 자기 예측에 미치는 영향이며, 이상치 진단에 씁니다. 정사영 행렬의 대각 성분이고 합이 입니다.
심화 6. 분해들을 한자리에 정리하고 언제 무엇을 쓸지 정하세요.
풀이. 지금까지 배운 분해와 앞으로 배울 것을 정리합니다.
| 분해 | 형태 | 조건 | 비용 | 강의 |
|---|---|---|---|---|
| LU | 정사각 | 70 | ||
| 콜레스키 | A=CC^ | 대칭 양정치 | 70 | |
| QR | 열 독립 | 2mn^ | 81 | |
| 고유분해 | A=S\Lambda S^ | 대각화 가능 | \sim n^ | 85 |
| 스펙트럼 | A=Q\Lambda Q^ | 대칭 | \sim n^ | 86 |
| SVD | A=U\Sigma V^ | 언제나 | \sim mn^ | 88 |
선택 기준을 정리합니다.
| 하려는 일 | 쓸 분해 |
|---|---|
| 정사각계 풀기 | LU |
| 대칭 양정치계 | 콜레스키 |
| 최소제곱 | QR |
| 계수 부족한 최소제곱 | SVD |
| 거듭제곱, 미분방정식 | 고유분해 |
| 대칭행렬의 구조 | 스펙트럼 |
| 저계수 근사 | SVD |
| 조건수 진단 | SVD |
아래로 갈수록 강력하고 비쌉니다. SVD가 언제나 되지만 상수가 크므로, 필요한 만큼만 씁니다.
공통 전략은 70강에서 세운 것입니다.
| 분해 | 한 번 | 여러 번 |
|---|---|---|
| LU | 마다 n^ | |
| QR | 2mn^ | 마다 |
| 고유분해 | n^ | 거듭제곱은 n^ |
| SVD | mn^ | 근사와 의사역행렬 |
분해가 이 과목의 실용적 결론입니다. 앞의 단원들이 개념을 세웠고, 05단원과 06단원이 그 개념을 계산 가능한 형태로 만듭니다.
이 강의에서는 numpy만 씁니다. 수정 그람슈미트로 QR을 직접 만들어 라이브러리와 맞추고, 최소제곱을 두 방법으로 풀어 안정성을 비교하며, QR 알고리즘이 고유값에 수렴함을 확인합니다.
import numpy as np
def mgs_qr(A):
"""수정 그람슈미트로 축소형 QR을 계산합니다."""
m, n = A.shape; V = A.astype(float).copy()
Q = np.zeros((m,n)); R = np.zeros((n,n))
for j in range(n):
R[j,j] = np.linalg.norm(V[:,j]); Q[:,j] = V[:,j] / R[j,j]
for k in range(j+1, n):
R[j,k] = Q[:,j] @ V[:,k]; V[:,k] -= R[j,k] * Q[:,j]
return Q, R
# --- 문제 1: A = QR -----------------------------------------------------
A = np.column_stack([[1.,1.,0.],[1.,0.,1.],[0.,1.,1.]])
Q, R = mgs_qr(A)
print(np.round(Q,8).tolist())
# [[0.70710678, 0.40824829, -0.57735027], [0.70710678, -0.40824829, 0.57735027],
# [0.0, 0.81649658, 0.57735027]]
print(np.round(R,8).tolist())
# [[1.41421356, 0.70710678, 0.70710678], [0.0, 1.22474487, 0.40824829],
# [0.0, 0.0, 1.15470054]]
print(bool(np.allclose(Q@R, A)), bool(np.allclose(Q.T@Q, np.eye(3)))) # True True
print(bool(np.allclose(R, np.triu(R))), np.round(np.diag(R),8).tolist())
# True [1.41421356, 1.22474487, 1.15470054]
# --- 문제 1: R = Q^T A 이고 det 는 대각의 곱 ----------------------------
print((np.round(Q.T@A, 8)+0.0).tolist()) # R 과 같습니다
print("%.8f %.8f" % (abs(np.linalg.det(A)), abs(np.prod(np.diag(R)))))
# 2.00000000 2.00000000
# --- 문제 2: 축소형과 완전형 --------------------------------------------
B = np.column_stack([[1.,1.,1.,1.],[0.,1.,2.,3.]])
Qb, Rb = mgs_qr(B)
print(Qb.shape, Rb.shape, bool(np.allclose(Qb@Rb, B))) # (4, 2) (2, 2) True
print(np.round(Rb,8).tolist()) # [[2.0, 3.0], [0.0, 2.23606798]]
Qf, Rf = np.linalg.qr(B, mode='complete')
print(Qf.shape, Rf.shape, bool(np.allclose(Qf@Rf, B))) # (4, 4) (4, 2) True
print(bool(np.allclose(Qf.T@Qf, np.eye(4)))) # True
# --- 문제 3: QR 로 최소제곱 ---------------------------------------------
b = np.array([1., 3., 2., 5.])
xq = np.linalg.solve(Rb, Qb.T@b) # R x = Q^T b
xn = np.linalg.solve(B.T@B, B.T@b) # 정규방정식 (대조용)
print(np.round(xq,8).tolist(), np.round(xn,8).tolist()) # [1.1, 1.1] [1.1, 1.1]
r = b - B@xq
print((np.round(B.T@r,10)+0.0).tolist()) # [0.0, 0.0] (잔차가 좌영공간)
print("%.8f" % float(np.linalg.norm(r))) # 1.64316767
# --- 문제 5: 조건수를 제곱하는가 ----------------------------------------
for e in [1e-3, 1e-5, 1e-7]:
X = np.array([[1.,1.],[e,0.],[0.,e]])
xt = np.array([1., -1.]); bb = X @ xt
xn_ = np.linalg.solve(X.T@X, X.T@bb)
Qx, Rx = np.linalg.qr(X)
xq_ = np.linalg.solve(Rx, Qx.T@bb)
print("%.0e %.3e %.3e %.3e" % (e, np.linalg.cond(X),
np.linalg.norm(xn_-xt), np.linalg.norm(xq_-xt)))
# 1e-03 1.414e+03 1.014e-10 2.220e-16
# 1e-05 1.414e+05 1.171e-07 0.000e+00
# 1e-07 1.414e+07 1.131e-03 2.220e-16
# --- 심화 6: 비용 비교 (QR 대 정규방정식) -------------------------------
for (m, n) in [(1000,100), (1000,500), (5000,200)]:
print(m, n, "%.3e %.3e" % (2*m*n**2 - 2*n**3/3, m*n**2 + n**3/3))
# 1000 100 1.933e+07 1.033e+07
# 1000 500 4.167e+08 2.917e+08
# 5000 200 3.947e+08 2.027e+08
# --- 심화 2: QR 알고리즘이 고유값을 줍니다 ------------------------------
M = np.array([[4.,1.,0.],[1.,3.,1.],[0.,1.,2.]])
Ak = M.copy()
for _ in range(60):
Qk, Rk = np.linalg.qr(Ak); Ak = Rk @ Qk # 직교 닮음이라 고유값 보존
print(np.round(np.sort(np.diag(Ak))[::-1], 8).tolist())
print(np.round(np.sort(np.linalg.eigvalsh(M))[::-1], 8).tolist())
# [4.73205081, 3.0, 1.26794919]
# [4.73205081, 3.0, 1.26794919]
실행하면 주석과 같은 값이 나옵니다. 다섯 곳을 짚어 둡니다.
첫째, 직접 만든 이 를 만족하고 이 위삼각이며 대각이 모두 양수입니다. **80강의 그람슈미트 결과가 그대로 **입니다.
둘째, 가 확인되고 가 대각의 곱 와 일치합니다. 79강의 부피 배율이 여기서 다시 나옵니다.
셋째, 축소형은 와 이고 완전형은 와 입니다. **완전형에서만 **이며, 축소형 는 직교행렬이 아닙니다.
넷째, QR과 정규방정식이 같은 답 을 주고 잔차가 정확히 좌영공간에 있습니다. 이 예는 조건이 좋아 두 방법이 일치합니다.
다섯째가 이 강의의 핵심입니다. 조건이 나빠지면 갈립니다. 에서 정규방정식 오차가 이고 QR은 입니다. 이 배 작아질 때마다 정규방정식 오차가 약 배가 되어 법칙과 맞습니다.
심화 2의 QR 알고리즘도 인상적입니다. 분해와 곱을 번 반복했을 뿐인데 고유값이 소수 여덟 자리까지 정확합니다.
코드로 할 수 없는 일도 분명히 해 둡니다. QR 분해의 유일성은 예로 증명되지 않습니다. 심화 1의 논증이 그 자리를 맡습니다. 그리고 QR 알고리즘이 수렴하는 이유는 확인하지 않았습니다. 여기서는 수렴한다는 사실만 관측했고, 근거는 거듭제곱법과의 관계에 있으며 이 과목의 범위를 넘습니다. 또 이 예는 대칭행렬이라 잘 수렴했으며, 일반 행렬에서는 복소 고유값 때문에 실수 삼각형으로 수렴하지 않을 수 있습니다.
정답.
| 기호 | 읽는 법 | 뜻 |
|---|---|---|
| QR 분해 | 직교 곱하기 위삼각입니다 | |
| 축소형 | thin QR | 가 입니다 |
| 완전형 | full QR | 가 직교행렬입니다 |
| r_ | 대각 성분 | 독립성의 정도입니다 |
| 열 피벗팅 QR | 대각이 감소합니다 | |
| QR 알고리즘 | QR iteration | 고유값을 구합니다 |
| 기븐스 회전 | Givens | 성분 하나씩 없앱니다 |
| 레버리지 | leverage | 의 대각입니다 |
다음 82강에서는 최소제곱을 정면으로 다룹니다. 이 강의에서 QR로 푸는 방법을 보였는데, 그것이 왜 최소인지는 74강의 직교 관계에서 나옵니다. 정규방정식 를 유도하고, 계산에는 여전히 QR을 쓰는 이유를 정리합니다.