Trăn Gram-Schmidt

Xem các phiên từ Hội nghị chuyên đề WiML về các mô hình khuếch tán với KerasCV, ML trên thiết bị, v.v. Xem theo yêu cầu

• TenorFlow
• Tài nguyên
• xác suất
• API

Thực hiện thuật toán chỉnh hình Gram-Schmidt đã sửa đổi

tfp.math.gram_schmidt(
    vectors, num_vectors=None
)

Ở đây ta giả sử rằng các vectơ độc lập tuyến tính. Các vectơ không sẽ không thay đổi, nhưng cũng sẽ sử dụng một phép lặp đối với num_vectors

Từ 1]. "MGS tương đương về mặt số với phân tích thừa số QR của Chủ nhà được áp dụng cho ma trận A được tăng cường bằng ma trận vuông không có phần tử nào ở trên cùng. "

Ghi chú lịch sử, xem [1]. Gram-Schmidt "đã sửa đổi" được tạo ra bởi Laplace [2], để loại bỏ và không phải là một thuật toán trực giao. Gram-Schmidt "cổ điển" thực sự đến sau [2]. Gram-Schmidt cổ điển đôi khi có sự mất trực giao nghiêm trọng đối với các ma trận có điều kiện kém, điều này sẽ được thảo luận thêm trong [1]

Người giới thiệu

[1] Bjorck, A. (1994). Số trực giao gram-schmidt. Đại số tuyến tính và các ứng dụng của nó, 197, 297-316

[2] P. S. Laplace, Lý thuyết phân tích xác suất. Bổ sung hàng đầu, thưa bà. Tòa án, Paris, 1816

[3] E. Schmidt, über die Auflosung linearer Gleichungen mit unendlich vielen Unbekannten, Rend. lưu thông. Chiếu. Pulermo (1) 25. 53-77 (1908)

Việc thực hiện được lấy từ Smile. https. //github. com/haifengl/smile/blob/master/core/src/main/java/smile/projection/RandomProjection. java

double[][] X;
Math.unitize(X[0]);
for (int i = 1; i < p; i++) {
    for (int j = 0; j < i; j++) {
        double t = -Math.dot(X[i], X[j]);
        Math.axpy(t, X[j], X[i]);
    }
    Math.unitize(X[i]);
}

ở đâu

• Math.unitize unit-chuẩn hóa vector
• Math.axpy cập nhật một mảng bằng cách thêm bội số của một mảng khác y = a * x + y
• ở đây X[i] là cột thứ $i$ của $X$

con trăn

def normalize(v):
    return v / np.sqrt(v.dot(v))

n = len(A)

A[:, 0] = normalize(A[:, 0])

for i in range(1, n):
    Ai = A[:, i]
    for j in range(0, i):
        Aj = A[:, j]
        t = Ai.dot(Aj)
        Ai = Ai - t * Aj
    A[:, i] = normalize(Ai)

thừa số QR

Chúng ta có thể nghĩ về Quy trình Gram-Schmidt bằng ngôn ngữ ma trận (như đối với Loại bỏ Gaussian đưa chúng ta đến Nhân tố hóa LU)

Trong đại số tuyến tính, các cơ số trực giao có nhiều tính chất hay. Ví dụ, các ma trận bao gồm các vectơ cột trực giao (a. k. một. ma trận trực giao) có thể dễ dàng đảo ngược bằng cách chuyển đổi ma trận. Ngoài ra, chẳng hạn, việc chiếu các vectơ trên các không gian con được kéo dài bởi các vectơ trực giao với nhau sẽ dễ dàng hơn. Quá trình Gram-Schmidt là một thuật toán quan trọng cho phép chúng ta chuyển đổi một cơ sở tùy ý thành một cơ sở trực giao bao trùm cùng một không gian con. Trong bài đăng này, chúng tôi sẽ triển khai và trực quan hóa thuật toán này ở dạng 3D với thư viện Nguồn mở phổ biến manim

Quy trình Gram-Schmidt

Xét một vài không gian con ddd tùy ý AA=⟨a1​,…,ad​⟩

được kéo dài bởi các vectơ độc lập tuyến tính a1,…,ada_1, \dots, a_da1< . Chúng tôi muốn tìm cơ sở sau bao trùm cùng một không gian con. ​,…,ad​. We want to find the following basis spanning the same subspace: A=⟨q1,…,qd⟩A = \langle q_1, \dots, q_d \rangleA=⟨q1​,…,qd​⟩

Quy trình Gram-Schmidt là một thuật toán lập trình động điển hình, vì ý tưởng cốt lõi đằng sau nó là tạo ra ⟨q1,…,qi⟩\langle q_1, \dots, q_i . Công trình thi công như sau. ⟨q1​,…,qi​⟩ an orthogonal basis assuming ⟨q1,…,qi−1⟩\langle q_1, \dots, q_{i-1} \rangle⟨q1​,…,qi−1​⟩ is already such a basis. The construction works as follows:

Bước 1. Đối với q1q_1q1​ . a1a_1a1​ and normalize it: q1. =a1∥a1∥q_1. = \frac{a_1}{\lVert a_1 \rVert}q1​:=∥a1​∥a1​​

Bước 2. q2q_2q2​ phải . Chúng ta có thể làm điều này bằng cách trừ đi hình chiếu của ⟨q1⟩\langle q_1 \rangle⟨q1​⟩, so we need to find out the component of a2a_2a2​ that is orthogonal to ⟨q1⟩\langle q_1 \rangle⟨q1​⟩. We can do this by subtracting the projection of a2a_2a2​< . onto ⟨q1⟩\langle q_1 \rangle⟨q1​⟩ from a2a_2a2​: q2′. =a2−(q1⊤a2)q1q1⊤q1=a2−(q1⊤a2)q1q'_2. = a_2 - \frac{(q_1^\top a_2) q_1}{q_1^\top q_1} = a_2 - (q_1^\top a_2) q_1q2′​:=a2​−q1⊤​q1​(q1⊤​a2​)q1​​=a2​−(q1⊤​a2​)q1​

Dạng thứ hai đúng vì chúng tôi đã chuẩn hóa q1q_1q1​ vector and the denominator is thus equal to one.

Bước này có thể được hình dung như sau

Cuối cùng, chúng ta cần chuẩn hóa véc-tơ kết quả. q2. =q2′∥q2′∥q_2. = \frac{q'_2}{\lVert q'_2 \rVert}q2​:=∥q2′​∥q2′​​

Bước 3. q3q_3q3​ phải . ⟨q1,q2⟩\langle q_1, q_2 \rangle⟨q1​,q2​⟩, which is already an orthogonal basis, so q3q_3q3​ must be orthogonal to both ⟨q1⟩\langle q_1 \rangle⟨q1​⟩ and ⟨q2⟩\langle q_2 \rangle⟨q2​⟩ which means that we need to subtract the projections of a3a_3a3​ onto both ⟨q1⟩\langle q_1 \rangle⟨q1​⟩ and ⟨q2⟩\langle q_2 \rangle⟨q2​⟩ from a3a_3a3​: q3′. =a3−(q1⊤a3)q1−(q2⊤a2)q2q'_3. = a_3 - (q_1^\top a_3) q_1 - (q_2^\top a_2) q_2q3′ . ​:=a3​−(q1⊤​a3​)q1​−(q2⊤​a2​)q2​

Như trong bước trước, chúng ta chuẩn hóa véc-tơ. q3. =q3′∥q3′∥q_3. = \frac{q'_3}{\lVert q'_3 \rVert}q3​:=∥q3′​∥q3′​​

Tổng quát hơn, chúng ta có thể viết công thức cho bước 1≤i≤d1 \le i \le d1≤i≤d as follows:

Bước i. qiq_iqi​ phải . Bởi vì trong các bước trước chúng ta đã thực hiện ⟨q1,…,qi−1⟩\langle q_1, \dots, q_{i-1} \rangle⟨q1​,…,qi−1​⟩. Because in the previous steps we made ⟨q1,…,qi−1⟩\langle q_1, \dots, q_{i-1} \rangle⟨ . jq1​,…,qi−1​⟩ an orthogonal basis, it follows that qiq_iqi​ must be orthogonal to each vector from {qj:j{qj​:j<i} . Do đó, chúng tôi trừ các phép chiếu của aia_iai​ onto qjq_jqj​ for all jj<i from aia_iai​: qi′. =ai−(q1⊤ai)q1−(q2⊤ai)q2−⋯−(qi−1⊤ai)qi−1q'_i. = a_i - (q_1^\top a_i) q_1 - (q_2^\top a_i) q_2 - \dots - (q_{i-1}^\top a_i) q_{i-1}qi′​:=ai​−(q1⊤​ai​)q1​−(q2⊤​ai​)q2​−⋯−(qi−1⊤​ai​)qi−1​

Bây giờ chúng ta chỉ cần chuẩn hóa véc-tơ cơ sở mới. qi. =qi′∥qi′∥q_i. = \frac{q'_i}{\lVert q'_i \rVert}qi​:=∥qi′​∥qi′​​

Triển khai & Trực quan hóa

Chúng ta có thể triển khai thuật toán trực giao hóa Gram-Schmidt trong Python theo cách sau

import numpy as np


def gram_schmidt(A):
    
    (n, m) = A.shape
    
    for i in range(m):
        
        q = A[:, i] # i-th column of A
        
        for j in range(i):
            q = q - np.dot(A[:, j], A[:, i]) * A[:, j]
        
        if np.array_equal(q, np.zeros(q.shape)):
            raise np.linalg.LinAlgError("The column vectors are not linearly independent")
        
        # normalize q
        q = q / np.sqrt(np.dot(q, q))
        
        # write the vector back in the matrix
        A[:, i] = q

Các vectơ cơ sở đầu vào được cung cấp dưới dạng các cột của ma trận, sau đó được thuật toán sửa đổi tại chỗ. Điều quan trọng là phải xác minh rằng qi′q'_iqi′ . ​ is not zero, because if qi′=0q'_i = 0qi′​=0, then it means that ai∈⟨q1,…,qi−1⟩=⟨a1,…,ai−1⟩a_i \in \langle q_1, \dots, q_{i-1} \rangle = \langle a_1, \dots, a_{i-1} \rangleai​∈⟨q1​,…,qi−1​⟩=⟨a1​,…,ai−1​⟩ and ⟨a1,…,ad⟩\langle a_1, \dots, a_d \rangle⟨a1​,…,ad​⟩ is therefore not a basis (vectors are not linearly independent).

Để trực quan hóa thuật toán, tôi đã thêm một vài câu lệnh yield vào phần triển khai ở trên, vì chúng cho phép chúng tôi tách biệt thuật toán khỏi mã trực quan hóa một cách tao nhã. Để trực quan hóa, tôi sẽ sử dụng thư viện manim phổ biến

# file: "gram-schmidt.py"
from manimlib.imports import *
from enum import Enum


class Action(Enum):
    UPDATE_MATRIX_REMOVE_Q = 1
    ADD_PROJECTION = 2
    REMOVE_PROJECTIONS_SET_Q = 3
    NORMALIZE_Q = 4


def gram_schmidt(A):
    
    (n, m) = A.shape
    
    for i in range(m):
        
        q = A[:, i] # i-th column of A
        
        for j in range(i):
            yield (Action.ADD_PROJECTION, np.dot(A[:, j], A[:, i]) * A[:, j])
            q = q - np.dot(A[:, j], A[:, i]) * A[:, j]
        
        if np.array_equal(q, np.zeros(q.shape)):
            raise np.linalg.LinAlgError("The column vectors are not linearly independent")
        
        yield (Action.REMOVE_PROJECTIONS_SET_Q, q)
        
        # normalize q
        q = q / np.sqrt(np.dot(q, q))

        yield (Action.NORMALIZE_Q, q)
        
        # write the vector back in the matrix
        A[:, i] = q

        yield (Action.UPDATE_MATRIX_REMOVE_Q, None)


class GramSchmidt(ThreeDScene):

    CONFIG = {
        "x_axis_label": "$x$",
        "y_axis_label": "$y$",
        "basis_i_color": GREEN,
        "basis_j_color": RED,
        "basis_k_color": GOLD,
        "q_color": PURPLE,
        "q_shifted_color": PINK,
        "projection_color": BLUE
    }

    def create_matrix(self, np_matrix):

        m = Matrix(np_matrix)

        m.scale(0.5)
        m.set_column_colors(self.basis_i_color, self.basis_j_color, self.basis_k_color)

        m.to_corner(UP + LEFT)

        return m

    def construct(self):
        
        M = np.array([
            [1.0, 1.0, -3.0],
            [3.0, 2.0, 1.0],
            [-2.0, 0.5, 2.5]
        ])

        axes = ThreeDAxes()

        axes.set_color(GRAY)

        axes.add(axes.get_axis_labels())

        self.set_camera_orientation(phi=55 * DEGREES, theta=-45 * DEGREES)

        basis_helper = TextMobject("$ a_1 $", ",", "$ a_2 $", ",", "$ a_3 $")
        basis_helper[0].set_color(self.basis_i_color)
        basis_helper[2].set_color(self.basis_j_color)
        basis_helper[4].set_color(self.basis_k_color)

        q_helper = TextMobject("Current $ q_i $ vector")
        q_helper[0].set_color(self.q_color)

        q_shifted_helper = TextMobject("Shifted $ q_i $ vector")
        q_shifted_helper[0].set_color(self.q_shifted_color)

        projection_helper = TextMobject("Projection vectors")
        projection_helper[0].set_color(self.projection_color)

        helper = VGroup(
            projection_helper,
            q_helper,
            q_shifted_helper,
            basis_helper
        )

        helper.arrange(
            DOWN,
            aligned_edge = RIGHT,
            buff=0.1
        )

        helper.to_corner(UP + RIGHT)

        self.add_fixed_in_frame_mobjects(helper)

        # matrix

        matrix = self.create_matrix(M)

        self.add_fixed_in_frame_mobjects(matrix)

        # axes & camera

        self.add(axes)

        self.begin_ambient_camera_rotation(rate=0.15)

        i_vec = Vector(M[:, 0], color=self.basis_i_color)
        j_vec = Vector(M[:, 1], color=self.basis_j_color)
        k_vec = Vector(M[:, 2], color=self.basis_k_color)

        self.play(
            GrowArrow(i_vec),
            GrowArrow(j_vec),
            GrowArrow(k_vec),
            Write(helper)
        )

        self.wait()

        projection_vectors = []

        q = None

        for (action, payload) in gram_schmidt(M):

            if action == Action.UPDATE_MATRIX_REMOVE_Q:

                assert not q is None

                M_rounded = np.round(M.copy(), 2)
                
                self.remove(matrix)

                matrix = self.create_matrix(M_rounded)

                self.add_fixed_in_frame_mobjects(matrix)

                i_vec_new = Vector(M[:, 0], color=self.basis_i_color)
                j_vec_new = Vector(M[:, 1], color=self.basis_j_color)
                k_vec_new = Vector(M[:, 2], color=self.basis_k_color)

                animation_time = 2.0

                self.play(
                    FadeOut(q, run_time=animation_time * 0.75),
                    ReplacementTransform(i_vec, i_vec_new, run_time=animation_time),
                    ReplacementTransform(j_vec, j_vec_new, run_time=animation_time),
                    ReplacementTransform(k_vec, k_vec_new, run_time=animation_time)
                )

                self.wait()

                i_vec, j_vec, k_vec = i_vec_new, j_vec_new, k_vec_new

            elif action == Action.ADD_PROJECTION:

                p = Vector(payload, color=self.projection_color)

                projection_vectors.append(p)

                self.play(GrowArrow(p))

                self.wait()

                if len(projection_vectors) == 2:

                    first_projection_end = projection_vectors[0].get_end()

                    p_shifted = Arrow(first_projection_end, first_projection_end + payload, buff=0, color=self.projection_color)

                    projection_vectors[1] = p_shifted

                    self.play(ReplacementTransform(p, p_shifted))

                    self.wait()

            elif action == Action.REMOVE_PROJECTIONS_SET_Q:

                if not projection_vectors:

                    q = Vector(payload, color=self.q_color)

                    self.play(GrowArrow(q))

                    self.wait()

                    continue

                last_projection_end = projection_vectors[len(projection_vectors) - 1].get_end()

                q_shifted = Arrow(last_projection_end, last_projection_end + payload, buff=0, color=self.q_shifted_color)

                self.play(GrowArrow(q_shifted))

                self.wait()

                q = Vector(payload, color=self.q_color)

                self.play(
                    ReplacementTransform(q_shifted, q),
                    *[FadeOut(p) for p in projection_vectors]
                )

                self.wait()

                projection_vectors = []

            elif action == Action.NORMALIZE_Q:

                q_normalized = Vector(payload, color=self.q_color)

                self.play(ReplacementTransform(q, q_normalized))

                self.wait()

                q = q_normalized

            else:
                assert False

            self.wait(1)

        assert np.allclose(M.T @ M, np.identity(3))

        self.wait(15)

Bằng cách hiển thị hoạt hình với manim gram-schmidt.py GramSchmidt -m -p, chúng tôi nhận được. Rất tiếc, trình duyệt của bạn không hỗ trợ thẻ video

gam là gì

Tại sao là gam

Tại sao Gram biến đổi

Ai đã nghĩ ra Gram