Microsoft MVP성태의 닷넷 이야기
Math: 22. 행렬로 바라보는 피보나치 수열 [링크 복사], [링크+제목 복사]
조회: 4619
글쓴 사람
홈페이지
첨부 파일
 

행렬로 바라보는 피보나치 수열

다음의 책을 보니 재미있는 내용이 있습니다. ^^

프로그래머를 위한 선형대수
; http://www.yes24.com/24/goods/39446808

(평을 보시면 아시겠지만, 저 역시 추천하고 싶은 책입니다. ^^)

249페이지에 보면 "자기회귀모델(AR: AutoRegressive)"의 이산시간에 대한 예로,

오늘의 ζ(t)는 어제의 ζ(t - 1), 이틀 전의 ζ(t - 2), 사흘 전의 ζ(t - 3)과 오늘의 u(t)에 따라 다음과 같이 정해진다.

ζ(t) = -0.5ζ(t - 1) + 0.34ζ(t - 2) + 0.08ζ(t - 3) + 2u(t)

초기 조건 ζ(0) = 0.78, ζ(-1) = 0.8, ζ(-2) = 1.5

소개가 되면서 다음과 같이 행렬 표현을 합니다.




저걸 보니, 피보나치 수열이 생각났습니다.

황금비율 증명 - 피보나치 수와 연분수의 관계
; https://www.sysnet.pe.kr/2/0/1312

역시 초깃값이 주어지고 x(t)는 x(t - 1)에 의해 결정되니까요. 따라서 위와 같은 기준으로 피보나치 수열을 바라보면 다음과 같이 정리가 됩니다.

ζ(t) = 1ζ(t - 1) + 1ζ(t - 2)

초기 조건 ζ(0) = 1, ζ(-1) = 0

간단하게 t = 1 ~ 4까지 테스트하면 이렇게 되고,

t = 1일 때, ζ(1) = ζ(1 - 1) + ζ(1 - 2) = ζ(0) + ζ(-1) = 1 + 0 = 1
t = 2일 때, ζ(2) = ζ(2 - 1) + ζ(2 - 2) = ζ(1) + ζ(0) = 1 + 1 = 2
t = 3일 때, ζ(3) = ζ(3 - 1) + ζ(3 - 2) = ζ(2) + ζ(1) = 2 + 1 = 3
t = 4일 때, ζ(4) = ζ(4 - 1) + ζ(4 - 2) = ζ(3) + ζ(2) = 3 + 2 = 5

이를 행렬로 표현하면 다음과 같습니다.




마찬가지로 t = 1 ~ 4까지에 대해 행렬로 계산하면 이렇게 됩니다.







따라서 (초깃값 2개를 넘어) n 번째 피보나치 수열은,




보는 바와 같이 행렬 [1 1; 1 0]에 대해 n 승을 하고 그 값을 [1 0] 행렬에 곱하면 n 번째 피보나치 수열이 구해지는 것입니다. 실제로 octave 같은 도구를 이용해 다음과 같이 행렬 계산을 바로 해볼 수 있습니다.

function fib_1()
  
a = [1 1; 1 0]
b = [1;0]
a ^ 1 * b
a ^ 2 * b
a ^ 3 * b
a ^ 4 * b
a ^ 5 * b

endfunction

위의 함수를 실행하면 2*1 행렬이 5개가 출력되는 데 그것의 첫 번째 원소들을 보면 1, 2, 3, 5, 8로 피보나치 수열이 나옵니다.




행렬로 표현된 피보나치 계산에서 고윳값/고유벡터를 이용해 풀어보면 재미있는 결과가 나옵니다.

[선형대수학 #3] 고유값과 고유벡터 (eigenvalue & eigenvector)
; http://darkpgmr.tistory.com/105

행렬 [1 1; 1 0]에 대한 고윳값, 고유벡터를 계산해 보면,




위의 행렬식을 구하면,

= (1 - λ)(0 - λ) - 1
= λ2 - λ -1


위와 같이 구한 특성 다항식을 특성 방정식에 따라 0 값이 나오는 해를 구하면,

det(A - λ E) = 0

λ2 - λ -1 = 0

근의 공식에 따라,






와 같이 계산됩니다. 고윳값을 구했으니 고유벡터까지 구해볼까요? ^^



연립 방정식으로 풀으면,

(1 - λ)vx + vy = 0
vx - λvy = 0
vx = λvy

따라서, vx가 vy의 λ배로 이뤄진 무수히 많은 벡터 = [λt, t]


그럼 고유 벡터를 아무거나 다음과 같이 선정할 수 있습니다.



따라서 고윳값 λ의 2가지 값에 대해,






이 중에서 고유 벡터를 [(1 + sqrt(5)) / 2, 1]인 쌍으로 골라 보겠습니다. 이를 다시 Gram-Schmidt 정규 직교로 바꾸면,

Matlab/Octave로 Gram-Schmidt 정규 직교 집합 구하는 방법
; https://www.sysnet.pe.kr/2/0/11235

(0.52573, -0.85065), (-0.85065, -0.52573)로 구할 수 있습니다. 즉, 이 2개의 벡터 각각에 대응하는 λ배의 모든 벡터들이 고유 벡터들이 됩니다.




실제로 위의 과정들을 간단하게 octave로 구할 수 있습니다.

a = [1 1; 1 0]
[ev, ei] = eig(a)

ev = 
    0.52573 -0.85065
   -0.85065 -0.52573

ei =

Diagonal Matrix

   -0.61803   0
   0          1.61803

또한, Av = λv인 것도 다음과 같이 쉽게 계산해볼 수 있습니다.

a * [0.52573, -0.85065]'
ans =

  -0.32492
   0.52573

-0.61803 * [0.52573 -0.85065]'
ans =

  -0.32492
   0.52573




피보나치 수열의 고윳값과 고유벡터를 구했으니 n 번째 값을 구하는 방법에 대해 행렬의 성질로 다시 살펴보겠습니다.

"[선형대수학 #3] 고유값과 고유벡터 (eigenvalue & eigenvector)" 글에 보면 다음과 같은 공식이 나옵니다.

A = 행렬
P = 행렬 A의 고유벡터들을 열벡터로 하는 행렬
Λ = 교윳값들을 대각 원소로 하는 대각 행렬

AP = PΛ
A = PΛP-1

이를 기반으로 A의 n 승을 다음과 같이 쉽게 구할 수 있는 방법을 포함하고 있습니다.

Ak = (PΛP-1)k
   = (PΛP-1)(PΛP-1)......(PΛP-1)
   = PΛkP-1
   = Pdiag(λk1,......,λkn)P-1

따라서, 가령 5번째 피보나치 수를 구하고 싶다면 고유 벡터와 그것의 역행렬만 구한 후 고윳값 2개를 대각 행렬로 갖는 것만 5 승을 해주면 되는 것입니다. 이것을 octave로 다음과 같이 테스트할 수 있습니다.

ev * ei ^ 5 * inverse(ev)
ans =

    8.0000  5.0000
    5.0000  3.0000

즉, 고윳값을 알기 전에는 다음과 같은 행렬 계산이었지만,




고윳값을 알게 된 이상, 그것은 대각행렬의 n 승으로 바뀌었기 때문에 단순히 스칼라 값인 고윳값 2개만 n 승을 해주면 되는 문제로 바뀐 것입니다.




[이 글에 대해서 여러분들과 의견을 공유하고 싶습니다. 틀리거나 미흡한 부분 또는 의문 사항이 있으시면 언제든 댓글 남겨주십시오.]





[최초 등록일: ]
[최종 수정일: 9/11/2017 ]

Creative Commons License
이 저작물은 크리에이티브 커먼즈 코리아 저작자표시-비영리-변경금지 2.0 대한민국 라이센스에 따라 이용하실 수 있습니다.
by SeongTae Jeong, mailto:techsharer@outlook.com

비밀번호

댓글 쓴 사람
 




... 31  32  33  34  35  36  37  [38]  39  40  41  42  43  44  45  ...
NoWriterDateCnt.TitleFile(s)
11328정성태10/15/20176666.NET Framework: 692. C# - 하나의 바이너리로 환경에 맞게 32비트/64비트 EXE를 실행하는 방법파일 다운로드1
11327정성태10/15/20173732.NET Framework: 691. AssemblyName을 .csproj에서 바꾼 경우 빌드 오류 발생하는 문제파일 다운로드1
11326정성태10/15/20174156.NET Framework: 690. coreclr 소스코드로 알아보는 .NET 4.0의 모듈 로딩 함수 [1]
11325정성태10/14/20173584.NET Framework: 689. CLR 4.0 환경에서 DLL 모듈의 로드 주소(Base address) 알아내는 방법
11324정성태10/18/20174044디버깅 기술: 101. windbg - "*** WARNING: Unable to verify checksum for" 경고 없애는 방법
11322정성태10/18/20173375디버깅 기술: 100. windbg - .NET 4.0 응용 프로그램의 Main 메서드에 Breakpoint 걸기
11321정성태10/11/20173701.NET Framework: 688. NGen 모듈과 .NET Profiler
11320정성태10/12/20173531.NET Framework: 687. COR_PRF_USE_PROFILE_IMAGES 옵션과 NGen의 "profiler-enhanced images"
11319정성태10/11/20176159.NET Framework: 686. C# - string 배열을 담은 구조체를 직렬화하는 방법
11318정성태10/7/20174474VS.NET IDE: 122. 비주얼 스튜디오에서 관리자 권한을 요구하는 C# 콘솔 프로그램 제작
11317정성태10/4/20178870VC++: 120. std::copy 등의 함수 사용 시 _SCL_SECURE_NO_WARNINGS 에러 발생
11316정성태7/9/20204202디버깅 기술: 99. (닷넷) 프로세스(EXE)에 디버거가 연결되어 있는지 아는 방법 [3]
11315정성태9/29/201713529기타: 68. "시작하세요! C# 6.0 프로그래밍: 기본 문법부터 실전 예제까지" 구매하신 분들을 위한 C# 7.0/7.1 추가 문법 PDF [8]
11314정성태9/28/20173964디버깅 기술: 98. windbg - 덤프 파일로부터 닷넷 버전 확인하는 방법
11313정성태9/25/20173858디버깅 기술: 97. windbg - 메모리 덤프로부터 DateTime 형식의 값을 알아내는 방법파일 다운로드1
11312정성태9/25/20174734.NET Framework: 685. C# - 구조체(값 형식)의 필드를 리플렉션을 이용해 값을 바꾸는 방법파일 다운로드1
11311정성태9/20/20173376.NET Framework: 684. System.Diagnostics.Process 객체의 명시적인 해제 권장
11310정성태9/19/20174798.NET Framework: 683. WPF의 Window 객체를 생성했는데 GC 수집 대상이 안 되는 이유 [3]
11309정성태9/18/20183217개발 환경 구성: 335. Octave의 명령 창에서 실행한 결과를 복사하는 방법
11308정성태9/13/20173976VS.NET IDE: 121. 비주얼 스튜디오에서 일부 텍스트 파일을 무조건 메모장으로만 여는 문제파일 다운로드1
11307정성태9/13/20176669오류 유형: 421. System.Runtime.InteropServices.SEHException - 0x80004005
11306정성태10/4/20184487.NET Framework: 682. 아웃룩 사용자를 위한 중국어 스팸 필터 Add-in
11305정성태9/12/20175394개발 환경 구성: 334. 기존 프로젝트를 Visual Studio를 이용해 Github의 신규 생성된 repo에 올리는 방법 [1]
11304정성태9/11/20173638개발 환경 구성: 333. 3ds Max를 Hyper-V VM에서 실행하는 방법
11303정성태9/11/20174906개발 환경 구성: 332. Inno Setup 파일의 관리자 권한을 제거하는 방법
11302정성태9/11/20173440개발 환경 구성: 331. SQL Server Express를 위한 방화벽 설정
... 31  32  33  34  35  36  37  [38]  39  40  41  42  43  44  45  ...