실습 · 연습문제 2개

대칭 양의 정부호(SPD) 행렬 $\mathbf{A}$에 대한 선형시스템 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 conjugate gradient(CG)로 푼다. CG는 $\mathbf{A}$의 원소가 아니라 행렬–벡터 곱 $\mathbf{v} \mapsto \mathbf{A}\mathbf{v}$만 쓴다.

$\mathbf{x}_0 = \mathbf{0}$, $\mathbf{r}_0 = \mathbf{b}$, $\mathbf{p}_0 = \mathbf{r}_0$에서 시작해 다음을 반복한다.

$$\alpha_k = \frac{\mathbf{r}_k^\top\mathbf{r}_k}{\mathbf{p}_k^\top\mathbf{A}\mathbf{p}_k}, \quad \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k\mathbf{p}_k, \quad \mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k\mathbf{A}\mathbf{p}_k,$$ $$\beta_k = \frac{\mathbf{r}_{k+1}^\top\mathbf{r}_{k+1}}{\mathbf{r}_k^\top\mathbf{r}_k}, \quad \mathbf{p}_{k+1} = \mathbf{r}_{k+1} + \beta_k\mathbf{p}_k.$$

함수 conjugate_gradient(matvec, b, tol=1e-10, max_iter=1000)를 완성한다.

• matvec은 벡터 v를 받아 A @ v를 돌려주는 함수다. 행렬 자체는 주어지지 않는다.
• 상대 잔차 $\|\mathbf{r}_k\| / \|\mathbf{b}\|$가 tol 이하가 되면 멈춘다. 시작할 때 이미 만족하면 반복하지 않는다.
• (x, n_iter)를 반환한다. n_iter는 수행한 반복(위 갱신) 횟수다.
• matvec은 반복마다 한 번만 호출한다(테스트가 호출 횟수를 센다). $\mathbf{A}\mathbf{p}_k$를 변수에 담아 $\alpha_k$와 $\mathbf{r}_{k+1}$ 계산에 함께 쓴다.

정확한 산술에서 CG는 $n \times n$ 시스템을 최대 $n$번 안에 푼다. 테스트는 결과의 정확도, 반복 수, matvec 호출 수를 확인한다.

CSR(compressed sparse row) 형식은 $m \times n$ 희소 행렬을 세 배열로 저장한다.

• data: 0이 아닌 값들을 행 순서대로 이어 붙인 배열 (길이 nnz)
• indices: 각 값의 열 번호 (길이 nnz)
• indptr: 행 $i$의 값이 data[indptr[i]:indptr[i+1]]에 있음을 알려 주는 포인터 (길이 $m+1$)

예를 들어

$$\mathbf{A} = \begin{bmatrix} 4 & 0 & 1 \\ 0 & 0 & 0 \\ 2 & 3 & 0 \end{bmatrix}$$

은 data = [4, 1, 2, 3], indices = [0, 2, 0, 1], indptr = [0, 2, 2, 4]다. 둘째 행은 비어 있어 indptr[1] == indptr[2]이다.

함수 csr_matvec(data, indices, indptr, x)를 완성해 $\mathbf{y} = \mathbf{A}\mathbf{x}$를 반환한다.

• 행의 수는 len(indptr) - 1이다. 반환 shape은 (m,)이다.
• scipy.sparse나 dense 행렬을 만들지 않고 세 배열만으로 계산한다.
• 빈 행의 결과는 0이다.

힌트: 행 $i$에 대해 data[indptr[i]:indptr[i+1]]와 x[indices[indptr[i]:indptr[i+1]]]의 내적을 구한다. 계산량은 nnz에 비례하고 $m \times n$과는 무관하다.