rep movsd

Cholesky분해는 양정부호행렬을 하삼각행렬과 그 전치행렬의 곱으로 표현하는 것이다라는 내용은 기초 선형대수학 교과서에 나오는 얘기다. 수학과 담을 쌓고 사는 사람이 아니면 늦어도 학부 3학년에는 이것을 배운다. 수학적 중요성은 잠깐 잊어버리더라도 계산상으로 Cholesky분해는 꽤 중요하다. 우선 LU분해에 걸리는 시간의 절반밖에 안걸린다. 게다가, 수치상 안정적이어서 다른 알고리즘에 비해 계산과정상 누적되는 오차의 영향을 덜 받는다. 단점도 물론 있다. 자세한 것은 선형대수 교과서나 수치해석 교과서를 보면 잘 나온다.

주요한 사용 예를 들자면: 최적화 문제의 답을 수치해석적으로 구해야만 하는 경우에 목적함수가 터무니없이 이상하게 생기지 않았으면 가중최소제곱법을 반복해서 답을 구할 수 있다. 이때 가중최소제곱법의 답을 구할때 Cholesky분해를 이용한다. (물론, 다른 방법으로 가중최소제곱법의 답을 구할 수도 있다. 예를 들어 MINPACK은 QR분해를 이용한다.)

실제로 계산하는 법은 의외로 간단하다. 이미 구현되어 제공되는 함수들도 많다. 예를들어 LAPACK의 DPOTRF는 배정도실수 행렬을 다루는 Fortran서브루틴이다. 같은 기능을 하는 함수가 NAG와 IMSL에도 물론 있다. 그런데, 이론적으로 Choleksy분해는 양반정부호행렬까지 확장될 수 있다. 그러나, 대부분의 라이브러리의 Cholesky분해 함수들은 양정부호행렬만 다룬다. 양반정부호행렬을 다루는 구현이 이미 1968년에 공개되었다는 점을 고려하면 이런 상황이 좀 의아하다. 내가 아는한 IMSL의 imsl_f_lin_sol_nonnegdef()만이 양반정부호행렬을 다루는 C 구현이다.

이하의 구현은 위에 링크된 1968년에 공개된 Fortran 구현을 C로 다시 구현하면서 Farebrother and Berry의 버그수정과 그것을 고려한 Barrett and Healy의 개선을 포함한 것이다. Healy의 원래 구현은 입력행렬과 출력행렬이 같은 주소를 가리켜도 되었는데, Barrett and Healy의 개선점을 포함하면서 그것이 허용되지 않게 되었다. 양정부호행렬만 분해한다면 이 함수는 LAPACK의 DPPTRF와 같은 기능을 한다.

함수 원형

int cholesky(double A[], int n, double U[], int *nullity)

A (입력)

양반정부호 행렬. 하삼각행렬부분만 1차원배열로 저장한다. 즉 A는 n*(n+1)/2개의 원소를 A₁₁, A₂₁, A₂₂, A₃₁, ... 순으로 저장한 1차원 배열.

n (입력)

A의 차원. 즉, 수학적으로 A는 n×n행렬.

U (출력)

Cholesky분해의 결과로 나오는 하삼각행렬. 역시 1차원 배열로 저장된다.

nullity (출력)

A의 영공간의 차원. NULL이면 출력이 없다. NULL을 전달했으나 이후에 영공간의 차원을 알고자 한다면, 대각원소들 중에 0.0인 것들의 갯수를 세면 된다.

0

작업 성공

-1

비정상적인 매개변수: n <= 0이거나 A와 U가 같은 주소를 가리키는 경우.

1 ... n

A가 양반정부호행렬이 아니며, 리턴값에 해당하는 행에서 처음으로 양반정부호행렬이 아님이 발견됨

구현
수식자체와 양정부호행렬임을 검사하는 방법은 아무 선형대수 교과서나 수치해석 교과서를 봐도 잘 설명되어 있다. 위에 링크된 논문들과 그에 따른 구현에서 중요한 것은 언제 어떻게 선형종속임을 결정하는가 하는 부분이다. (약간 줄 수가 많아 접어놓음.) C 코드를 보려면 클릭

#include <math.h>
#define eta 1.0e-9 /* Healy (1968)에서 설정된 값 */

int
cholesky(double A[], int n, double U[], int *nullity)
{
    int     i, j, k, l, m, icol, irow;
    double  w, x;

    if (n <= 0 || A == U)
	return -1;

    if (nullity)
	*nullity = 0;

    icol = j = k = 0;
    do {
	x = eta * eta * A[k + icol];
	irow = 0, l = 0;
	do {
	    w = A[k];
	    for (i = irow, m = j; i > 0; --i)
		w -= U[l++] * U[m++];
	    if (irow == icol)
		break;
	    if (U[l] == 0.0) {
		if (w * w > fabs(x * A[l]))
		    return k + 1;
		U[k] = 0.0;
	    }
	    else
		U[k] = w / U[l];
	    ++l, ++k;
	} while (++irow <= icol);

	if (fabs(w) > fabs(eta * A[k])) {
	    if (w < 0.0)
		return k + 1;
	    U[k] = sqrt(w);
	}
	else {
	    U[k] = 0.0;
	    if (nullity)
		++(*nullity);
	}
	++k, ++icol;
	j += icol;
    } while (icol < n);
    return 0;
}

더 해볼만한 것
원 논문의 바로 다음에 이어지는 논문에는 이 함수를 이용해서 양(반)정부호행렬의 (일반화)역행렬을 구하는 함수가 공개되어 있다. 같은 기능을 하는 것을 스스로 만들어쓰려면 BLAS의 DTPSV를 이용하면 된다.

실행파일 제공하는 곳은 cygwin, mingw, gnuwin32 등등 여러군데 있는데, 하나같이 불필요한 DLL을 설치하기를 강요하거나, 패치파일이나 패치대상 파일의 개행문자가 자기들이 정한대로 따르지 않으면 그냥 자빠져 버리는 등, 사용성에 있어서 참 짜증스럽게 만들어 놓았다.

어차피 patch는 소스가 공개되어 있으니 좀 고치면 원하는대로 바꿀 수 있다. 목표는, 우선 개행문자의 제약을 없애고, 다음으로 유닉스에서처럼 --posix옵션이 주어지지 않으면 백업파일을 자동으로 만들게 하는 것이다. 어차피 patch를 이용해 패치하는 파일들은 소스코드들이니, 개행문자를 가지고 패치가 엉터리라고 까다롭게 구는 것은 별 의미가 없다고 생각된다.

이 패치를 제작하는 과정에서 컴파일러는 마이크로소프트에서 무료로 제공하는 것을 썼다. 곁다리: gcc를 쓸거면 어차피 cygwin이나 mingw을 설치해야 gcc를 쓸 수 있으니, 거기서 제공하는 patch실행파일을 쓰는게 여기의 것을 쓰는 것보다 일관성의 측면에서 낫다. 게다가 mingw용 patch는 reject파일도 unified diff로 만드는 옵션이 있다. context diff에 익숙하지 않으면 꽤 유용한 옵션이다.

소스는 ftp://ftp.gnu.org/gnu/patch/patch-2.5.9.tar.gz에서 받아오면 되고, 패치파일은 아래에서 받으면 된다.

명령행 요약: (패치파일 맨처음과 같은 내용. 현재 스킨에서는 마지막 줄의 일부만 보인다.)

tar xzf patch-2.5.9.tar.gz
cd patch-2.5.9
patch -p0 -b < patch.patch
cl /I. /DHAVE_CONFIG_H /Ded_PROGRAM=\"ed\" /Fepatch.exe error.c addext.c argmatch.c backupfile.c basename.c dirname.c getopt.c getopt1.c inp.c maketime.c partime.c patch.c pch.c quote.c quotearg.c quotesys.c util.c version.c xmalloc.c

---

@ 패치를 패치한다고 하니까 말이 꼬인다.

invalid-file

MMX 지원이 없는 고전 펜티엄에서 빠른 memcpy()를 구현하기 위해 다음과 같은 코드를 쓰던 적이 있었다.

L0:     fild    QWORD PTR [esi]
        fistp   QWORD PTR [edi]
        add     esi,8
        add     edi,8
        sub     ecx,8
        ja      L0

그런데, 정수를 FPU 레지스터로 읽어들이는 fild는 부동소수점 실수를 읽어들이는 fld보다 느리다. 메모리로 쓰는 fistp 역시 fstp보다 느리다. 고전 펜티엄에서는 그 차이가 두드러진다. 어셈블리 초보의 입장에서 보면, 64비트 정수나 배정도 실수나 모두 64비트 크기를 가지는 자료이니, 조금 더 빠른 fld/fstp의 조합을 쓰는 것이 당연하게 보인다. 그런 생각을 실제로 벤치마크로 보여주는 코드도 웹에 떠돌아 다닌다.

링크된 벤치마크 코드를 돌려봐도 잘만 돌아가는데, 왜 이것이 초보나 저지르는 실수인가? 그것은, 코드자체에 있는 것이 아니라, 부동소수점 자료의 표현때문이다. fld/fstp를 쓰는 memcpy() 구현은 보통때에는 별 이상없이 작동하지만, 복사하는 내용이 특정한 값이면 프로그램 실행이 멎어버리는 경우가 있다. 인텔 매뉴얼 1권의 용어로는, SNaN 값이 복사대상이면 그런 일이 벌어진다. (엄밀하게 말하자면, 실행이 멈출지 아닐지는 프로그램에서 예외처리를 어떻게 하는가와 운영체제의 예외처리 정책에 따른다. 프로그램에서 FPU 예외처리를 하지 않고 운영체제의 기본 예외처리를 따른다면, 운영체제의 정책에 따라 실행이 멎을 수도 있고 계속 실행할 수도 있다. 프로그램에서 FPU 예외처리를 하면, 실행은 계속 할 수 있겠지만, NaN 값을 FPU 레지스터로 읽어들이는 경우에는 보통 실수 값을 읽을 때보다 한참 느리게 작동한다. 그러므로, 예외발생/예외처리/NaN처리의 과정을 거치면 엄청난 속도상의 손실이 발생한다.)

이것을 직접 확인해 볼 수 있는 예로서, 위의 링크에 있는 파일을 받아서 begin() 함수를 다음과 같이 약간 고쳐서 실행해 볼 수 있다. (빨간 글씨로 되어 있는 부분을 추가한다.)

void begin( LPBYTE mem1, LPBYTE mem2, int size, const char* text ) {
	memset(mem1,0x55,size);
	memset(mem2,0xff,size); // 0xAA를 0xff로 바꾸고 다음 두 줄을 추가.
	for (int i = 0; i < size; i += 8)
	    mem2[i+6] = 0xf0;
	printf(text);
}

인텔 32비트 CPU에서 조건부 분기 없이 32비트 정수의 절대값 구하기는 이미 잘 알려진 코드다. 32비트 운영체제가 대세가 되면서 한동안 어셈블리 코딩을 하지 않던 나에게 다시 인텔 어셈블리어 코딩을 하도록 끌어들인 그 코드. 조건부 분기하는 것보다 빠른 것은 물론 더 짧다.

; 입력:  eax = 절대값을 구할 정수
; 출력:  eax = 입력값의 절대값
; 기타 사용 레지스터:  edx
    cdq
    xor eax,edx
    sub eax,edx

초보를 위한 설명:

• CDQ명령으로 EAX 레지스터의 부호를 EDX로 복사하면, EDX는 -1이거나 0이다.
• XOR연산을 0에 대해 하면, 바뀌는 것 없다. 2의 보수로 표현한 -1에 대해 하면 NOT연산과 같다.
• SUB명령에 인수로 0을 주면, 바뀌는 것 없다. -1을 주면 INC명령과 같다.

어느정도 인텔 어셈블리 코딩을 해온 사람이라면 이 코드를 하도 많이 써서 관용구로 외우고 있을 것이다. 매크로로 만들어 놓은 사람도 분명 있을 것이고... 게다가, 앞뒤로 한줄씩만 더하면 바로 C표준함수 labs()이다. 물론, 이 코드를 내가 본 것이 벌써 10년전이니까, 요즘 컴파일러들은 labs()를 이렇게 바로 인라인할 것이라고 추측된다.

누군가 어셈블리어를 검색해서 찾아온 것을 보고 나서 32비트 초보 시절의 기억을 되새겨 보았다.

• 2의 보수를 사용하는 기계에서, 정수 i의 맨 아래쪽 1 비트(1로 세트된 비트 중 LSB에 제일 가까운 비트)를 0으로 만들려면 C 수식으로 i & (i-1)를 쓴다. 부호없는 정수로 생각하면, 1을 빼는 연산이 항상 맨 아래 비트를 0으로 만들고, 그보다 더 LSB쪽에 있는 비트는 1로 만든다는 것을 쉽게 알 수 있다.

  이 방법을 사용한 예는 Ken Thompson이 Unix에 썼다고 하는 전설이 있는 비트 수 세기 루프이다.

  /* i에 비트 수를 세고자 하는 값이 들어있다. */
  bitcount = 0;
  while (i != 0) {
      bitcount++;
      i &= i - 1;
  }

  이 코드는 2의 보수를 사용하는 모든 기계에서 작동한다.

• 반대로, 2의 보수를 사용하는 기계에서 정수 i의 맨 아래쪽 1 비트만 남기고 다 지우려면, C 수식으로 i & -i를 쓴다. 이것도 2의 보수방식으로 부호를 바꾸는 연산이 NOT연산후에 1을 증가시킨다는 것을 기억하면 이 코드를 이해할 수 있다.

  이 방법을 응용하면, 어떤 정수가 2의 누승인지 쉽게 알아낼 수 있다. 이진법에서 2의 누승은 단 하나의 비트만 1이고 나머지 비트는 모두 0이다. 그러므로, 이렇게 하면 빨리 알아낼 수 있다:

  if ((i & -i) == i) { /* 2의 누승 */ }

  이 역시 머신워드 크기에 관계없이 2의 보수를 사용하는 모든 기계에서 작동한다.

• 맨 아래쪽 1 비트가 몇번째 비트인지 알아내는 것은 루프를 돌아도 알아낼 수 있지만, 인텔계열 386이상의 CPU에서는 bsf명령을 쓰면 한번에 알 수 있다. BSD계통의 Unix에는 같은 목적으로 ffs()라는 함수가 제공된다.