Assignment 8

6 minute read

과제 8

과제 8
점수 0
제출 False

선형 공간상의 임의의 linearly independent 벡터들이 주어졌다고 해보겠습니다. 우리는 Gram-Schmidt algorithm (GS)를 이용하여 주어진 벡터들로 span되는 선형 공간의 orthonormal basis를 얻어낼 수 있습니다. GS는 매우 간단하고 강력하지만, 실제로 구현할때는 실제 값과의 오차가 크게 나타날 수 있습니다. 심지어는 (정확하게) orthogonal한 벡터들을 얻어내지 못 할 수도 있습니다. 여기서는 이런점을 개선한 modified Gram-Schmidt algorithm을 소개하려 합니다.

1. Classical Gram-Schmidt algorithm (CGS)

선형 공간 \(V\)상에 선형 독립인 \(n\)개의 벡터 \(\{v_1,\,v_2,\,\cdots,\,v_n\}\)들이 주어졌다고 하겠습니다. 이 때,

\[\begin{align*} v'_1 &= v_1,\\ v'_2 &= v_2 - \mbox{proj}_{v_1'}v_2,\\ v'_3 &= v_3 - \mbox{proj}_{v_1'}v_3 - \mbox{proj}_{v_2'}v_3,\\ &\qquad\qquad \vdots\\ v'_n &= v_n - \sum_{i=1}^{n-1}\mbox{proj}_{v_i'}v_n,\\ u_i &= \frac{v_i'}{\lVert v'_i \rVert},\qquad i = 1, 2, \cdots, n \end{align*}\]

와같은 과정을 통해 얻어낸 \(\{u_1,\,u_2,\,\cdots,\,u_n\}\)들은 \(\mbox{span}\{v_1,\,v_2,\,\cdots,\,v_n\}\)의 orthonormal basis 입니다. 위 과정을 잘 살펴보면, \(v_k'\)를 생성할 때, \(v_k\)를 \(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\)로 span되는 평면상으로 projection한 뒤에 \(v_k\)에 빼는 과정으로 볼수 있습니다. 만약 ‘\(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\)로 span되는 평면’을 구할때부터 오차가 생긴다면, 우리가 구해낸 \(v_k'\)는 실제로 구해져야 하는 것과는 다른 방향을 가리킬 것 입니다. 다음 classical GS의 코드를 보겠습니다.

%--- The following is the function file 'CGramSchmidt.m'. ---%

% Find an orthonormal basis for col(A) when A has full column rank.
function Q = CGramSchmidt(A)

[m, n] = size(A);

% Initialize the matrix Q as an m*n zero matrix.
Q = zeros(m, n);
for i = 1 : n
    % normalization
    v = A(:, i);
    for j = 1 : i-1
        % Subtract each component of orthogonal projection of v
        % onto the subspace spanned by the vector Q(:, j).
        q = Q(:, j);
        v = v - (q' * A(:, i)) * q;
    end
    Q(:, i) = v / norm(v);
end
end
% Q is an m*n matrix whose columns form an orthonormal basis for col(A).

여기서 특히 오차가 많이 발생하는 부분은

v = v - (q' * A(:, i)) * q;

입니다. 즉, 각 \(v_i\)로 projection할 때 오차가 발생하게 되며 이 오차들은 쌓이고 쌓여서 ‘\(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\)로 span되는 평면’이 제대로 계산되지 않는 문제가 생기게 됩니다.

2. Modified Gram-Schmidt algorithm (MGS)

GS과정에서 생기는 오차를 줄이는 방법은 간단합니다. 각 projection을 단계적으로 시행하여, 이 전 projection에서 생기는 오차들을 다음 projection과정에서 삭제하는 것이죠. 다만, 단계적인 projection들을 모두 거친 결과물은 classical GS와 똑같아야 합니다.

\[\begin{align*} v'_{k, 1} &= v_k - \mbox{proj}_{u_1}v_k,\\ v'_{k, 2} &= v'_{k, 1} - \mbox{proj}_{u_2}v'_{k, 1},\\ v'_{k, 3} &= v'_{k, 2} - \mbox{proj}_{u_3}v'_{k, 2},\\ &\qquad\qquad \vdots\\ v'_{k, k-1} &= v'_{k, k-2} - \mbox{proj}_{u_{k-1}}v'_{k, k-2}\\ u_k &= \frac{v_{k, k-1}'}{\lVert v'_{k, k-1} \rVert}, \qquad k = 1,\,2,\cdots,\,n \end{align*}\]

위의 여러줄의 수식은 모두 하나의 \(u_k\)를 생성하는 과정입니다. Classical GS와는 달리, \(v_1\)으로 projection하여 얻은 결과를 다시 \(v_2\)로 projection하고, 또 다시 \(v_3\)로 반복하여 projection하고 있습니다.

연습문제

  1. MGS의 알고리즘을 참고하여, modified Gram-Schmidt를 수행하는 코드를 작성해보세요. (힌트: Exercise 7.13)

  2. Classical GS와 modified GS를 비교하세요.

    임의 행렬

    1. 임의의 \(100 \times 100\) 행렬 \(A\)를 생성하세요. rank명령어를 사용하여 \(A\)가 non-singular인 것을 확인하세요.
      (높은 확률로 \(A\)는 non-singular matrix입니다, 만약 \(A\)가 singular matrix라면 다시 생성하세요.)

       A = rand(100);
       fprintf('Rank of A : %d\n', rank(A));
      
    2. \(\mbox{col}(A)\)의 orthonormal basis를 구성하세요. CGS를 이용한 결과를 \(Q_1\), MGS를 이용한 결과를 \(Q_2\)라고 할 때, 단위행렬과의 오차를 비교하세요.

       E1 = norm(Q1' * Q1 - eye(100));
       E2 = norm(Q2' * Q2 - eye(100));
       fprintf('Error of CGS = %g\n', E1)
       fprintf('Error of MGS = %g\n', E2)
      

    안 좋은 행렬 (Ill-conditioned matrix)

    1. 수치적으로 좋지않은 행렬 \(B\)를 생성하세요:

       % Ex)
       B = 10^-5 * rand() * eye(100) + hilb(100);
      
    2. CGS를 이용해 얻은 결과를 \(P_1\), MGS를 이용한 결과를 \(P_2\)라고 할 때, 단위행렬과의 오차를 비교하세요.

       E1 = norm(P1' * P1 - eye(100));
       E2 = norm(P2' * P2 - eye(100));
       fprintf('Error of CGS = %g\n', E1)
       fprintf('Error of MGS = %g\n', E2)
      

Assignment 8

Assignment 8
Score 0
Submission False

Let’s say we are given random linearly independent vectors in a linear space. We can use the Gram-Schmidt algorithm (GS) to obtain the orthonormal basis of the linear space spanned by given vectors. GS is very simple and powerful, but when implemented, there may be a large error from the actual value. We may not even be able to get (exactly) orthogonal vectors. Here, we will show a modified Gram-Schmidt algorithm that improves this point.

1. Classical Gram-Schmidt algorithm (CGS)

Suppose that \(n\) linearly independent vectors \(\{v_1,\,v_2,\,\cdots,\,v_n\}\) are given in the linear space \(V\). We obtained \(\{u_1,\,u_2,\,\cdots,\,u_n\}\) through the following process:

\[\begin{align*} v'_1 &= v_1,\\ v'_2 &= v_2 - \mbox{proj}_{v_1'}v_2,\\ v'_3 &= v_3 - \mbox{proj}_{v_1'}v_3 - \mbox{proj}_{v_2'}v_3,\\ &\qquad\qquad \vdots\\ v'_n &= v_n - \sum_{i=1}^{n-1}\mbox{proj}_{v_i'}v_n,\\ u_i &= \frac{v_i'}{\lVert v'_i \rVert},\qquad i = 1, 2, \cdots, n \end{align*}\]

which is an orthonormal basis of \(\mbox{span}\{v_1,\,v_2,\,\cdots,\,v_n\}\). Looking closely at the above process, when creating \(v_k'\), \(v_k\) is projected to the plane spanned by \(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\) and subtracted to it self. If there is an error from finding ‘the plane spanned by \(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\)’, the \(v_k'\) that we found will point in a different direction than what should be. Let’s look at the classical GS code below.

%--- The following is the function file 'CGramSchmidt.m'. ---%

% Find an orthonormal basis for col(A) when A has full column rank.
function Q = CGramSchmidt(A)

[m, n] = size(A);

% Initialize the matrix Q as an m*n zero matrix.
Q = zeros(m, n);
for i = 1 : n
    % normalization
    v = A(:, i);
    for j = 1 : i-1
        % Subtract each component of orthogonal projection of v
        % onto the subspace spanned by the vector Q(:, j).
        q = Q(:, j);
        v = v - (q' * A(:, i)) * q;
    end
    Q(:, i) = v / norm(v);
end
end
% Q is an m*n matrix whose columns form an orthonormal basis for col(A).

Here, particularly where the error occurs is:

v = v - (q' * A(:, i)) * q;

In other words, an error occurs when projecting each \(v_i\), and these errors are accumulated. Thus there is a problem that the ‘plane spanned by \(\{v_1,\,v_2,\,\cdots,\,v_{k-1}\}\)’ is not calculated properly.

2. Modified Gram-Schmidt algorithm (MGS)

It is easy to reduce errors in the GS process. Project in stepwise. Then the errors that occur in the previous projection are removed in the next projection process. Simultaniously, the result of going through all of the step-by-step projections should be the same as the classical GS.

\[\begin{align*} v'_{k, 1} &= v_k - \mbox{proj}_{u_1}v_k,\\ v'_{k, 2} &= v'_{k, 1} - \mbox{proj}_{u_2}v'_{k, 1},\\ v'_{k, 3} &= v'_{k, 2} - \mbox{proj}_{u_3}v'_{k, 2},\\ &\qquad\qquad \vdots\\ v'_{k, k-1} &= v'_{k, k-2} - \mbox{proj}_{u_{k-1}}v'_{k, k-2}\\ u_k &= \frac{v_{k, k-1}'}{\lVert v'_{k, k-1} \rVert}, \qquad k = 1,\,2,\cdots,\,n \end{align*}\]

All of the above multi-line formulas are the process of creating a single \(u_k\). Unlike the Classical GS, the result obtained by projecting with \(v_1\) is projected to \(v_2\), and projected repeatedly with \(v_3\).

Exercises

  1. Write a code that executes modified Gram-Schmidt by referring to the algorithm of MGS. (Hint: Exercise 7.13)

  2. Compare Classical GS and Modified GS.

    Random matrix

    1. Generate a random \(100 \times 100\) matrix \(A\). Check that \(A\) is non-singular using the rank command.
      (There is a high probability that \(A\) is a non-singular matrix, if \(A\) is a singular matrix, recreate it.)

       A = rand(100);
       fprintf('Rank of A : %d\n', rank(A));
      
    2. Generate an orthonormal basis of \(\mbox{col}(A)\). Let the result using CGS be \(Q_1\) and the result using MGS be \(Q_2\). Compare the error with the unit matrix.

       E1 = norm(Q1' * Q1 - eye(100));
       E2 = norm(Q2' * Q2 - eye(100));
       fprintf('Error of CGS = %g\n', E1)
       fprintf('Error of MGS = %g\n', E2)
      

    Bad matrix (Ill-conditioned matrix)

    1. Create a matrix \(B\) that is not numerically good:

       % Ex)
       B = 10^-5 * rand() * eye(100) + hilb(100);
      
    2. When the result obtained using CGS is \(P_1\) and the result using MGS is \(P_2\), compare the error with the unit matrix.

       E1 = norm(P1' * P1 - eye(100));
       E2 = norm(P2' * P2 - eye(100));
       fprintf('Error of CGS = %g\n', E1)
       fprintf('Error of MGS = %g\n', E2)
      

Leave a comment