|
- # !/usr/bin/env python
- # -*- coding: utf-8 -*-
- """
- @Time : 2021/5/8 15:09
- @Author : Albert Darren
- @Contact : 2563491540@qq.com
- @File : RGB_split.py
- @Version : Version 1.0.0
- @Description : TODO
- @Created By : PyCharm
- """
- import numpy as np
-
- # 自己原创平方根法
- def cholesky_decomposition(sym_pos_define_matrix: np.ndarray, right_hand_side_vector: np.ndarray):
- rows, columns = sym_pos_define_matrix.shape
- # judge if it is a square matrix
- if rows == columns:
- # 判断是否对称
- if (sym_pos_define_matrix.T == sym_pos_define_matrix).all():
- for k in range(rows): # 判断各阶顺序主子式是否为0,即判定是否正定
- if np.linalg.det(sym_pos_define_matrix[:k + 1, :k + 1]) <= 0:
- raise Exception("cannot decompose")
- else:
- # 初始化单位下三角矩阵L
- square_ones_matrix = np.ones((rows, columns))
- lower_triangle_matrix = np.tril(square_ones_matrix)
- for j in range(rows):
- prod = 0
- for k in range(j):
- prod += lower_triangle_matrix[j, k] * lower_triangle_matrix[j, k]
- lower_triangle_matrix[j, j] = np.sqrt(sym_pos_define_matrix[j, j] - prod)
-
- for i in range(j + 1, rows):
- prod = 0
- for k in range(j):
- prod += lower_triangle_matrix[i, k] * lower_triangle_matrix[j, k]
- lower_triangle_matrix[i, j] = (sym_pos_define_matrix[i, j] - prod) / lower_triangle_matrix[j, j]
-
- else:
- raise Exception("error,the input matrix must be a symmetric-matrix")
- else:
- raise Exception("ERROR:please pass a square matrix.")
- return np.linalg.inv(lower_triangle_matrix.transpose()) @ np.linalg.inv(
- lower_triangle_matrix) @ right_hand_side_vector
-
-
- if __name__ == '__main__':
- import cProfile
- symmetric_positive_define_matrix2 = np.array([[4, 2, -4, 0, 2, 4, 0, 0],
- [2, 2, -1, -2, 1, 3, 2, 0],
- [-4, -1, 14, 1, -8, -3, 5, 6],
- [0, -2, 1, 6, -1, -4, -3, 3],
- [2, 1, -8, -1, 22, 4, -10, -3],
- [4, 3, -3, -4, 4, 11, 1, -4],
- [0, 2, 5, -3, -10, 1, 14, 2],
- [0, 0, 6, 3, -3, -4, 2, 19]], dtype=np.float64)
- column_vector = np.array([0, -6, 20, 23, 9, -22, -15, 45], dtype=np.float64).reshape((8, 1))
- # print(cholesky_decomposition(symmetric_positive_define_matrix2, column_vector))
- cProfile.run("cholesky_decomposition(symmetric_positive_define_matrix2, column_vector)",
- filename="result.out", sort="cumulative")
|