init

Author: lonelyprince7

.vscode/launch.json

@@ -0,0 +1,15 @@
+{
+    // 使用 IntelliSense 了解相关属性。 
+    // 悬停以查看现有属性的描述。
+    // 欲了解更多信息，请访问: https://go.microsoft.com/fwlink/?linkid=830387
+    "version": "0.2.0",
+    "configurations": [
+        {
+            "name": "Python: 当前文件",
+            "type": "python",
+            "request": "launch",
+            "program": "${file}",
+            "console": "integratedTerminal"
+        }
+    ]
+}

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 HongYu Zhang
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

README.md

@@ -0,0 +1,3 @@
+# 《数值分析》
+
+用Python+Numpy+Pytorch实现《数值分析》(Timothy Sauer著)的所有算法

chaper0.多项式求值/霍纳法多项式求值.py

@@ -0,0 +1,26 @@
+'''
+Descripttion: horn法则多项式求值
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-07 11:40:20
+LastEditors: ZhangHongYu
+LastEditTime: 2021-03-07 11:59:29
+'''
+import numpy as np
+
+def poly(order, coeff, x, basis):
+    if basis is None:
+        basis = np.zeros((order,), dtype=np.float32)
+    # 初始递推值
+    y = coeff[-1]
+    for i in reversed(range(order)):
+        y = y * (x - basis[i]) + coeff[i]
+    return y
+
+if __name__ == '__main__':
+    # 令所有基点为0的普通多项式
+    res1 = poly(4, np.array([-1, 5,-3, 3, 2]), 1/2, np.array([0, 0, 0, 0]))
+    print(res1)
+    # 基点不为0的三阶插值多项式
+    res2 = poly(3, np.array([1, 1/2, 1/2, -1/2]), 1, [0, 2, 3])
+    print(res2)

chapter1.求解方程根/不动点迭代.py

@@ -0,0 +1,19 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-08 15:13:40
+LastEditors: ZhangHongYu
+LastEditTime: 2021-03-08 15:21:11
+'''
+import numpy as np
+import math
+def fpi(g, x0 ,k): #迭代k次,包括x0在内共k+1个数
+    x = np.zeros(k+1,)
+    x[0] = x0
+    for i in range(1, k+1):
+        x[i] = g(x[i-1])
+    return x[k]
+if __name__ == '__main__':
+    res = fpi(lambda x: math.cos(x), 0, 10)
+    print(res)

chapter1.求解方程根/二分法.py

@@ -0,0 +1,25 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-08 14:48:29
+LastEditors: ZhangHongYu
+LastEditTime: 2021-03-08 15:01:59
+'''
+import numpy as np
+import math
+def binary(f, a, b, tol):
+    if f(a)*f(b) >=0 :
+        raise ValueError("f(a)*f(b)<0 not satisfied!")
+    while (b-a)/2 > tol:
+        c=(a+b)/2 # 即使a和b是int,此处c自动转float了
+        if f(c) == 0: #c是一个解，完成
+            break
+        if f(a)*f(c)<0 :
+            b = c
+        else:
+            a = c
+    return (a+b)/2
+if __name__ == '__main__':
+    res = binary(lambda x: x**3+x-1, 0, 1, 5*pow(10, -5))
+    print(res)

chapter1.求解方程根/尝试求梯度.py

@@ -0,0 +1,21 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-08 16:41:50
+LastEditors: ZhangHongYu
+LastEditTime: 2021-06-06 20:51:56
+'''
+import torch
+#这里我们对w,b求梯度，事实上只要requires_grad=True,我们也可以对x,y求梯度
+x = torch.ones(5)  # input tensor
+y = torch.zeros(3)  # expected output  
+w = torch.randn(5, 3, requires_grad=True)
+b = torch.randn(3, requires_grad=True)
+z = torch.matmul(x, w)+b
+# print(z)
+loss = torch.nn.functional.binary_cross_entropy_with_logits(z, y)
+print(loss)  # loss是一个标量tensor(3.2646)，只有这样才可反向传播
+loss.backward()
+print(w.grad,'\n', b.grad)
+print(w.grad!=None)

chapter1.求解方程根/尝试迭代2.py

@@ -0,0 +1,22 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-06-06 21:48:48
+LastEditors: ZhangHongYu
+LastEditTime: 2021-06-06 21:58:20
+'''
+import torch
+a = torch.tensor([2., 3.], requires_grad=True)
+b = torch.tensor([6., 4.], requires_grad=True)
+Q = 3*a**3 - b**2
+external_grad = torch.tensor([1., 1.])
+Q.backward(gradient=external_grad)
+a.grad.zero_()
+b.grad.zero_()
+with torch.no_grad():
+    a += 1
+Q = 5*a*b
+Q.backward(gradient=external_grad)
+print(a.grad)
+print(b.grad)

chapter1.求解方程根/牛顿法.py

@@ -0,0 +1,37 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-08 15:57:56
+LastEditors: ZhangHongYu
+LastEditTime: 2021-06-06 22:17:31
+'''
+import numpy as np
+import math
+import torch
+#x.grad为Dy/dx(假设Dy为最后一个节点)
+def newton(x0, k, f): #迭代k次,包括x0在内共k+1个数
+    # 初始化计算图参数
+    x = torch.tensor([x0], requires_grad=True)
+    for i in range(1, k+1):
+        # 每次迭代都要调用函数前向传播一次并重新求梯度，类似于网络训练，然后更新变量(这里是x)
+        # 前向传播，注意x要用新的对象，否则后面y.backgrad后会释放
+        y = f(x)
+        # 因为调用backward的对象自身不算梯度，相当于d y /dy = 1->sy.grad其实还是None?
+        # 或者令loss=y.sum() loss.backward()
+        # y.backward(gradient = torch.tensor([1.0]))
+        y.backward() # y.grad是None
+        # 更新参数
+        # 注意此处必须调用x.add原地修改，否则会生成一份copy，调用.grad得None
+        # 修改式需关闭求导，否则会对leaf变量就地修改未定义行为(修改是计算图的一部分?)
+        with torch.no_grad(): 
+            x.sub_(torch.divide(y, x.grad))   
+        x.grad.zero_() # 清空梯度，使下一轮建立新的计算图，否则因为backward释放资源下一轮再backward出错
+        #注意x.grad不能是0，否则要出错使g(x)/x.grad变为none
+        # 这里g(x)为cos(x)就不行，这样在x=0 时, -sin0=0，故x.grad = 0
+    return x.detach().numpy()[0]
+if __name__ == '__main__':
+    f = lambda x: x**3 + x - 1
+    x0 = 1.0
+    res = newton(x0, 10, f)
+    print(res)

chapter12特征值和奇异值/计算SVD.py

@@ -0,0 +1,14 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-05-22 21:43:41
+LastEditors: ZhangHongYu
+LastEditTime: 2021-05-22 21:45:24
+'''
+#暴力计算ATA,然后求AT和A的特征值
+
+
+
+
+#计算B=[0 AT;A 0]的特征值和特征向量

chapter2.方程组/2.1.高斯消元/朴素高斯消元.py

@@ -0,0 +1,59 @@
+'''
+Descripttion: 朴素法无法解决零主元和高斯消元问题，高斯消元都考虑方非奇异阵(肯定是方的)
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-03-08 17:26:11
+LastEditors: ZhangHongYu
+LastEditTime: 2021-05-29 16:22:24
+'''
+import numpy as np
+eps = 1e-6
+# 消去步骤
+def gaussion_elimination(A, b): #假设A是方阵，A.shape[0] == A.shape[1], python中变量是传拷贝，数组等对象是传引用
+    assert(A.shape[0] == A.shape[1])
+    for j in range(A.shape[1]): #消去第j列的数
+        # abs(A[j ,j])为要消去的主元
+        if abs(A[j, j]) < eps:
+            raise ValueError("zero pivot encountered!")  #无法解决零主元问题
+            return
+        # 消去主对角线以下的元素A[i, j]
+        for i in range(j+1, A.shape[0]):
+            mult_coeff = A[i, j]/A[j, j]
+            # 对这A中这一行都进行更新
+            for k in range(j, A.shape[1]):
+                A[i, k] = A[i, k] - mult_coeff * A[j, k]
+            b[i] = b[i] - mult_coeff * b[j] #二维的b取单dim1的索引即[1]这种含单个元素的列表
+
+def gaussion_putback(A, b):
+    x = np.zeros((A.shape[0], 1))
+    for i in reversed(range(A.shape[0])): #算出第i个未知数
+        for j in range(i+1, A.shape[1]):
+            b[i] = b[i] - A[i, j] * x[j]
+        x[i] = b[i] / A[i, i]
+    return x
+
+if __name__ == '__main__':
+    A = np.array(
+        [
+            [1, 2, -1],
+            [2, 1, -2],
+            [-3, 1, 1]
+        ]
+    )
+    b = np.array(
+        [
+            [3],
+            [3],
+            [-6]
+        ]
+    )
+    gaussion_elimination(A, b)
+    x = gaussion_putback(A, b)
+    print(x)
+    #print(A, "\n", b)
+    
+
+            
+
+            
+

chapter2.方程组/2.2.LU分解/LU分解及回代.py

@@ -0,0 +1,92 @@
+'''
+Descripttion: 
+Version: 1.0
+Author: ZhangHongYu
+Date: 2021-05-22 21:53:43
+LastEditors: ZhangHongYu
+LastEditTime: 2021-05-29 21:01:44
+'''
+import numpy as np
+from copy import deepcopy
+eps = 1e-6
+# 消去步骤
+def LU_decomposition(A): #假设A是方阵，A.shape[0] == A.shape[1], python中变量是传拷贝，数组等对象是传引用
+    assert(A.shape[0] == A.shape[1])
+    U = deepcopy(A)
+    L = np.zeros(A.shape, dtype=np.float32)
+    for j in range(U.shape[1]): #消去第j列的数
+        # abs(U[j ,j])为要消去的主元
+        if abs(U[j, j]) < eps:
+            raise ValueError("zero pivot encountered!")  #无法解决零主元问题
+            return
+        L[j, j] = 1
+        # 消去主对角线以下的元素A[i, j]
+        for i in range(j+1, U.shape[0]):
+            mult_coeff = U[i, j]/U[j, j]
+            L[i, j] = mult_coeff
+            # 对这A中这一行都进行更新
+            for k in range(j, U.shape[1]):
+                U[i, k] = U[i, k] - mult_coeff * U[j, k]
+            
+    return L, U
+
+#常规的上三角进行回代(此例中对角线不为0)
+def gaussion_putback_U(A, b):
+    x = np.zeros((A.shape[0], 1))
+    for i in reversed(range(A.shape[0])): #算出第i个未知数
+        for j in range(i+1, A.shape[1]):
+            b[i] = b[i] - A[i, j] * x[j]
+        x[i] = b[i] / A[i, i]
+    return x
+
+#下三角进行回代(此例中对角线不为0)
+def gaussion_putback_L(A, b):
+    x = np.zeros((A.shape[0], 1))
+    for i in range(A.shape[0]): #算出第i个未知数
+        for j in range(i):
+            b[i] = b[i] - A[i, j] * x[j]
+            #草,如果b矩阵初始化时是整形，3-6.99999976 = ceil(-3.99999) = -3，
+            # 直接给我向上取整(截断)约成整数了
+            # if i == A.shape[0] - 1:
+            #     print(A[i, j], "----", x[j], "----", A[i, j]*x[j])
+            #     print(b[i])
+        x[i] = b[i] / A[i, i]
+    return x
+
+def LU_putback(L, U, b):
+    # Ax = b => LUx = b ，令Ux = c
+    # 解 Lc = b
+    c = gaussion_putback_L(L, b) #上三角回代
+    print(c)
+    # 再解 Ux = c
+    x = gaussion_putback_U(U, c) #下三角回代
+    return x
+
+if __name__ == '__main__':
+    A = np.array(
+        [
+            [1, 2, -1],
+            [2, 1, -2],
+            [-3, 1, 1]
+        ],
+        dtype=np.float32
+    )
+    b = np.array(
+        [
+            [3],
+            [3],
+            [-6]
+        ],
+        dtype=np.float32 #注意，此处必须是浮点型，否则整形的话后面就自动舍入了
+    )
+    # 单纯的LU分解过程不会对b有影响
+    # 即消元与回代分离
+    
+    # 分解步骤
+    L, U = LU_decomposition(A) # A=LU
+    print(L)
+    print(U)
+    # 回代步骤
+    x = LU_putback(L, U, b)
+    print(x)
+    #print(A, "\n", b)