(App.1) Linear Solver

1. 개요

Linear Solver라고 하면,
통상 $Ax=b$ 형태의 수식에서 $x$를 구하는 방법을 말합니다.
앞으로 할 내용에서 기반을 이루는 부분이 되므로, 여기서는 간단한 방식에 대해서만 언급을 하고 넘어가도록 하겠습니다.

위 수식에서 $A$와 $b$는 행렬의 형태가 될 수 있으며,
$A$가 $(m\times n)$일 때, $b$는 $(m\times 1)$, $x$는 $(n\times 1)$ 입니다.

아래에 사용한 코드들의 경우는, 일반적으로 행렬 연산에 사용하는 Eigen라이브러리와 부가적으로 NVIDIA GPU를 사용하는 cuda기반의 cusolver를 이용해서 설명을 했습니다.

2. Solve (Least Square)

일반적으로 위의 수식에서 $A$와 $b$가 행렬의 경우는,
$A$의 역형렬($A^{-1}$)을 구해서, 양쪽에 곱해주면($A^{-1}Ax = A^{-1} b$), 바로 $x$가($x = A^{-1} b$) 나옵니다.

만약, $A$가 정방행렬이 아니라면($m\neq n$)라면 $A$의 Transpose$A^T$를 곱하여, 정방형렬로 만든 다음에 동일하게 역행렬을 구해서 곱하면 됩니다. 아래 처럼요.

$Ax=b \rightarrow A^TAx=A^Tb \rightarrow (A^TA)^{-1}A^TAx=(A^TA)^{-1}A^Tb \rightarrow x=(A^TA)^{-1}A^Tb$

이 방법은 잘 알려진 Least Square Method(최소자승법)입니다.
위에서 $A^TA$에 의해 각 항목들은 제곱 처리가 되기 때문에 Square라고 합니다.

$2\times 2$나 $3\times 3$의 경우와 같이 간단한 경우는 위의 방식으로 풀어도 무방합니다. 그런데 여기서 행렬의 크기가 커진다면, 역행렬을 구할 때, 꽤 많은 컴퓨팅 파워를 필요로 합니다. 즉 느립니다.

이를 해결하기 위해서, Linear Solver 교제들에서는 다양한 Solve방식을 제시합니다.

3. Linear Solver

Least Square방식의 경우, 역형렬을 구한 다음, 이를 곱해서 x를 구하게 됩니다. 그런데 역행렬을 구하는 과정에서 이루어지는 분해(decomposition)을 이용해서, 역행렬까지 구하지 않고, 바로 $A$를 단위행렬로 만드는 동시에, $b$행렬에도 동일한 연산을 수행하여 $x$를 구할 수 있습니다.

자세한 내용은 선형대수학 교재 등을 참고하시면 됩니다만, 굳이 정독하지 않아도 위에서 말한 Eigen과 같은 라이브러리에서 대부분을 다 지원하기 때문에 개념 정도만 알고 넘어가면 됩니다.

다시 돌아와서, decomposition을 보면, 다양한 방식이 존재합니다.
이는 $A$행렬의 형태에 따라서 아래와 같이 크게 3개로 분기가 됩니다.

먼저, Cholesky 의 경우는 정방행렬($m= n$)이면서도, 행렬이 Symmetry인 경우에만 적용할 수 있으며, 제한된 환경에서만 동작하는 만큼 제일 빠릅니다.

다행인 건, 최적화 문제의 경우 풀어야 되는 Matrix의 형태가 대부분 Symmetry이기에 이 방식을 주로 사용할 것입니다. 그래도 한번은 확인해야 겠죠..

두 번째는 LU방식으로 이는 $A$ 정방(Square)행렬이면서, Symmetry 가 아닌 경우입니다.

마지막, QR은 아예 행과 열 크기가 다른 행렬에 대한 경우입니다.

다시 정리하면 아래와 같습니다.

개별적으로 또 분기가 생기지만, 일단 우리가 풀고자하는 $A$의 형태에 따라서 이렇게 분기가 생긴다 정도만 알아두면 됩니다.

Solver	A 크기	Symmetry	속도
Cholesky decomposition	$m= n$	O	빠름
LU decomposition	$m= n$	X	보통
QR decomposition	$m\neq n$	X	느림

3. 소스

한번 코드를 작성해서 테스트를 해보겠습니다.

3.1. Eigen

Eigen::MatrixXd AMat, bMat;

//Load Data
...
...

//cholesky decomposition
Eigen::MatrixXd retChol = AMat.ldlt().solve(bMat);

//lu decomposition
Eigen::MatrixXd retLudc = AMat.fullPivLu().solve(bMat);

//qr decomposition
Eigen::MatrixXd retqrdc = AMat.fullPivHouseholderQr().solve(bMat);

//least square
Eigen::MatrixXd retLsts = AMat.inverse() * bMat;

3.2. 결과

아래는 $A$의 크기가 $1000\times 1000$인 경우에 대한 테스트 결과입니다.

결과를 보시면 제일 하단 Least Square로 계산것 보다, Cholesky로 푼 경우가 훨씬 빠른 것을 볼 수 있습니다. 만약 행렬의 크기가 2나 3정도이면, 수식이 바로 나오기 때문에 decomposition을 하지 않고 바로 계산이 가능하며 이때는 더욱 빨라집니다.

상세한 내용은 아래 Eigen문서 참조하세요.
https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html

Solver	시간(ms)
Cholesky decomposition	20
LU decomposition	228
QR decomposition	259
Least square	559

※ 참고로 Eigen의 경우, 비정방행렬A를 QR로 푸는 것 보다, Transpose를 통해 $A^TA$로 바꾸어 Cholesky로 푸는 것이 더 빠릅니다.

3.3 cusolver

CPU환경에서는 Eigen을 사용하면 되나, 행/열의 크기가 1000 이상을 넘어갈 경우, GPU를 사용하는 cusolver를 사용하게 되면 훨씬 더 빠른 결과를 얻을 수 있습니다. 채굴, 머신러닝이 아니라 이럴때도 GPU를 씁시다.

//parameter 선언
cusolverDnHandle_t handle = nullptr;
cudaStream_t       stream = nullptr;

cusolverStatus_t status = CUSOLVER_STATUS_SUCCESS;
cudaError_t      state;

//solver 사용을 위한 handle 및 stream생성
//첫 1회에 한해, 속도가 느림. 그 이후에는 내부적으로 만들어진 handle을 그냥 뱉어주는 것으로 보임.
status = cusolverDnCreate(&handle);
state = cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);
status = cusolverDnSetStream(handle, stream);

const int m = cols;
const int lda = m;
const int ldb = m;

double *Acopy = nullptr;
double *b = nullptr;

int     bufferSize = 0;
int *   info = NULL;
double *buffer = NULL;
int *   ipiv = NULL;

//cuda 메모리 할당
state = cudaMalloc((void **)&Acopy, sizeof(double) * m * m);
state = cudaMalloc((void **)&b, sizeof(double) * m);
state = cudaMalloc(&ipiv, sizeof(int) * m);

//cuda로 메모리 A, b 행렬 데이터 복사
state = cudaMemcpy(Acopy, aeMat.data(), sizeof(double) * m * m, cudaMemcpyHostToDevice);
state = cudaMemcpy(b, deMat.data(), sizeof(double) * m, cudaMemcpyHostToDevice);
cudaDeviceSynchronize();

//Cholesky Decomposition을 위한 buffer 크기 계산
cusolverDnDpotrf_bufferSize(handle, CUBLAS_FILL_MODE_LOWER, m, (double *)Acopy, lda,
                            &bufferSize);
cudaDeviceSynchronize();

//일부 파라미터 메모리 할당
cudaMalloc(&info, sizeof(int));
cudaMalloc(&buffer, sizeof(double) * bufferSize);
cudaMemset(info, 0, sizeof(int));
cudaDeviceSynchronize();

//Cholesky Decomposition
cusolverDnDpotrf(handle, CUBLAS_FILL_MODE_LOWER, m, Acopy, lda, buffer, bufferSize, info);
cusolverDnDpotrs(handle, CUBLAS_FILL_MODE_LOWER, m, 1, Acopy, lda, b, m, info);
cudaDeviceSynchronize();

//LU Decomposition
//cusolverDnDgetrf(handle, m, m, Acopy, lda, buffer, ipiv, info);
//cusolverDnDgetrs(handle, CUBLAS_OP_N, m, 1, Acopy, lda, ipiv, b, m, info);
//cudaDeviceSynchronize();

//계산 결과를 CPU로 다시 복사
state = cudaMemcpy(X0, b, sizeof(double) * m, cudaMemcpyDeviceToHost);

if (info) cudaFree(info);
if (buffer) cudaFree(buffer);
if (Acopy) cudaFree(Acopy);
if (b) cudaFree(b);
if (ipiv) cudaFree(ipiv);

동일 데이터로 테스트 시에
속도 향상이 보이며, 이는 행렬의 크기가 커질 수록 더욱 차이가 납니다. 그럼 얼마까지 빨라질까요?
사실 이 이후로는 속도보다는 메모리를 봐야 합니다. 빠르냐의 문제보다는 가능하냐의 문제로 보통 진행됩니다.
예를 들어, 50000x50000 행렬이 있다고 한다면, float으로 연산을 한다고 하더라도 50000x50000x4(바이트) = 약 9GB 가 필요합니다. 저같이 데이터를 하드코어하게 다루시는 분들은 이런 문제에 많이 직면들 하실겁니다. Eigen이라면 PC의 메모리를 cusolver라면 GPU의 메모리를 봐야 겠죠.

Solver	시간(ms)
Cholesky decomposition	15.110
LU decomposition	27.661

4. Sparse

앞서, 메모리 문제를 언급했습니다. 이를 언급한 이유는 바로 Sparse Matrix(희소행렬)를 설명하기 위해서입니다.
행렬 내부에서 요소가 대부분 0인 행렬을 Sparse Matrix라고 합니다.
대체로 위와 같이 큰 행렬들이 이러한 경향을 보입니다.

예를 들어 1000명의 사람들이 있을 때, 서로 알고 지내는 사람들끼리 짝을 행렬이 표현한다고 할 때, 1000x1000행렬에서 알고 지내는 사람들은 1로 두고, 모르면 0이라고 하겠습니다.
이 때, 사실 대부분 서로 잘 모를 가능성이 크기 때문에 행렬의 요소는 대부분 0이 됩니다. 이렇게 뭔가 관계를 나타내려고 할 때, 요소 간 관계가 없는 데이터가 많기 때문에 희소행렬이 자주 나타납니다.

이런 데이터는 필연적으로 위와 같은 메모리 낭비를 가져오기 때문에,
이를 해결하기 위해서 Sparse Matrix를 처리하기 위한 다양한 압축 기술들이 나와 있습니다. 가장 유명한 것이 CSR(Compressed sparse row)기법인데, 간단히 설명하면, 0이 아닌 데이터들만 모으고, 그 데이터의 행렬 상 위치를 찾아갈 수 있도록 하는 기법입니다.
자세한 건, 아래를 참고하세요.

https://ko.wikipedia.org/wiki/%ED%9D%AC%EC%86%8C%ED%96%89%EB%A0%AC

이런 행렬 구조의 장점은 메모리가 줄어드는 부분도 있으나,
또, 행렬 연산을 할 때, 유의미한 아이템에 대해서만 수행을 하기 때문에, 속도가 빨라질 수 있습니다.
그러나, 필연적으로 실제 행렬의 위치를 가져오기 위한 Conversion과정이 들어오게 되고, 이는 추가적인 부하로 작용하게 됩니다. 결론은 trade off가 발생한다는 겁니다. 행렬 요소 중에서 0의 비율이 어느 정도이냐에 따라, 일반적인 Dense matrix로 풀지, 아니면 Sparse로 풀지 결정하면 되는데, 저같은 경우는 경험적으로 0의 비율이 전체에서 7~80%이상인 경우에만 Sparse를 쓰는 것이 좋다고 보는 편입니다. 물론 정답은 아닙니다.

아래는 Eigen으로 Sparse Matrix를 정의해서 Linear Solver를 수행한 예시입니다.

Eigen::SparseMatrix<double>                    spaeMat = aeMat.sparseView();
Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(spaeMat);
Eigen::VectorXd                                        x = chol.solve(deMat.transpose());

위와 동일한 데이터로 Cholesky Decomposition을 이용해서 테스트 시에 성능을 보면 굉장히 빨라진 것을 볼 수 있는데, Eigen이 Sparse Matrix에 대해, 내부적으로 어떤 구조를 쓰는지는 모르겠으나, 뭔가 굉장히 최적화된 것으로 보입니다.

Type	시간(ms)
Dense	42.277
Sparse	5.414

아래는 CSR기법으로 Sparse Matrix를 만든 다음에
cusolver에 적용하여 Cholesky Decomposition을 이용해서 푼 소스입니다.

int m = aeMat.cols();

cusolverSpHandle_t handle = nullptr;
cusparseHandle_t   cusparseHandle = nullptr; 
cudaStream_t       stream = nullptr;
cusparseMatDescr_t descrA = nullptr;

int rowsA = m;           
int colsA = m;           
int nnzA = aeSMat.size();

int *   d_csrRowPtrA = nullptr;
int *   d_csrColIndA = nullptr;
double *d_csrValA = nullptr;
double *d_x = nullptr;
double *d_b = nullptr;

cusolverSpCreate(&handle);
cusparseCreate(&cusparseHandle);

cudaStreamCreate(&stream);
cusolverSpSetStream(handle, stream);

cusparseSetStream(cusparseHandle, stream);

cusparseCreateMatDescr(&descrA);

in00 = std::chrono::high_resolution_clock::now();

cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ZERO);

cudaMalloc((void **)&d_csrRowPtrA, sizeof(int) * (rowsA + 1));
cudaMalloc((void **)&d_csrColIndA, sizeof(int) * nnzA);
cudaMalloc((void **)&d_csrValA, sizeof(double) * nnzA);
cudaMalloc((void **)&d_x, sizeof(double) * colsA);
cudaMalloc((void **)&d_b, sizeof(double) * rowsA);

cudaMemcpy(d_csrRowPtrA, rowMat.data(), sizeof(int) * (rowsA + 1), cudaMemcpyHostToDevice);
cudaMemcpy(d_csrColIndA, colMat.data(), sizeof(int) * nnzA, cudaMemcpyHostToDevice);
cudaMemcpy(d_csrValA, aeSMat.data(), sizeof(double) * nnzA, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, deMat.data(), sizeof(double) * rowsA, cudaMemcpyHostToDevice);

double tol = 1.e-12;
int    reorder = 0;
int    singularity = 0;

//Cholesky Decomposition
cusolverSpDcsrlsvchol(handle, rowsA, nnzA, descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, d_b,
                        tol, reorder, d_x, &singularity);

성능을 보면, 미세하게 Dense가 더 빠른 것을 볼 수 있습니다. 이는 절대적이라기 보다는 Sparse한 정도에 따라 가변적입니다.
(놀라운 건, 이 정도 크기 데이터에서는 Eigen이 더 빠르네요.)

Type	시간(ms)
Dense	15.110
Sparse	16.728

5. 마침

이러저런 설명을 했지만, Cholesky와 LU, QR을 어떨 때 사용한다 정도만 숙지하면 됩니다. 나머지는 라이브러리가 다 알아서 해줍니다.