Levenberg-Macquardt python 代码

最新推荐文章于 2023-06-14 10:30:46 发布

Lingjie Fan

最新推荐文章于 2023-06-14 10:30:46 发布

阅读量554

点赞数 4

分类专栏： python 文章标签： python

本文链接：https://blog.csdn.net/weixin_40267373/article/details/109197300

版权

python 专栏收录该内容

2 篇文章 0 订阅

订阅专栏

Levenberg-Macquardt python 代码

介绍
代码
理解

介绍

作为一个拟合常用的算法，对于我来说可以用来调参使得模拟程序生成的光谱尽可能的匹配上实验测得的光谱， L-M算法网上能找到的 python 代码有好多bug，没办法只能自己写一个啦。

网上有bug的代码

先放代码，再聊聊对 Levenberg-Macquardt 算法、The Steepest Descent Method 以及 Newton Method 的理解。

代码

把L-M算法写成了一个类，相当于一个工具箱以后想用的时候直接import就行了。

import numpy as np


class Levenberg_Macquardt:
    Experiment_Data = None  # y
    Experiment_Wavelength = None # x
    Floating_Parameter = None
    Function = None
    k_max = 200

    def __init__(self):
        self.Experiment_Data = None
        self.Experiment_Wavelength = None
        self.Floating_Parameter = None  # Parameters must be float-type numbers.
        self.k_max = 200
        self.Function = None

    def Derivation(self, Wavelength, n):
        X1 = np.array(self.Floating_Parameter.copy(), dtype=np.float64)
        X2 = np.array(self.Floating_Parameter.copy(), dtype=np.float64)
        X1[n][0] -= 1e-6 * X1[n][0]
        X2[n][0] += 1e-6 * X2[n][0]
        P1 = self.Function(X1, Wavelength)
        P2 = self.Function(X2, Wavelength)
        D = (P2 - P1) / (2e-6 * np.array(self.Floating_Parameter, dtype=np.float64)[n][0])
        return D

    def Jaccobi(self):
        num_data = self.Experiment_Data.shape[1]
        num_params = self.Floating_Parameter.shape[0]
        J = np.mat(np.zeros((num_data, num_params)))
        for i in range(num_data):
            for j in range(num_params):
                J[i, j] = self.Derivation(self.Experiment_Wavelength[i], j)
        return J

    def Residual(self, x):
        num_data = self.Experiment_Data.shape[1]
        Residual = np.mat(np.zeros((num_data, 1)))
        for i in range(num_data):
            Residual[i] = self.Function(np.array(x), self.Experiment_Wavelength[i]) - self.Experiment_Data[0, i]
        return Residual

    def Converge(self):
        k, v = 0, 2
        x = self.Floating_Parameter
        tau = 1e-3
        epsilon1 = 1e-8
        epsilon2 = 1e-8

        J = self.Jaccobi()
        A = J.T * J
        f = self.Residual(x)
        g = J.T * f

        found = np.linalg.norm(g) <= epsilon1
        u = tau * A.max()

        while (not found) and k < self.k_max:
            H = A + u * np.eye(x.shape[0])
            h = -H.I * g

            if np.linalg.norm(h) <= epsilon2 * (np.linalg.norm(x) + epsilon2):
                found = True
            else:
                x_new = x + h

                f = self.Residual(x)
                f_new = self.Residual(x_new)
                F = 0.5 * np.linalg.norm(f) ** 2
                F_new = 0.5 * np.linalg.norm(f_new) ** 2
                gain_ratio = (F - F_new) / ((0.5 * h.T * (u * h - g))[0, 0])

                if gain_ratio > 0:
                    x = x_new
                    J = self.Jaccobi()
                    A = J.T * J
                    f = self.Residual(x)
                    g = J.T * f

                    found = np.linalg.norm(g) <= epsilon1

                    u = u * max([1 / 3, 1 - (2 * gain_ratio - 1) ** 3])
                    v = 2
                else:
                    u = u * v
                    v = 2 * v
        self.Floating_Parameter = x

下面是一个简单的测试程序，验证一下我的代码是可行的。

from Levenberg_Macquardt import Levenberg_Macquardt
import numpy as np

n = 100
wl = np.linspace(0,10,n)
a_true,b_true,c_true,d_true = -3,-10,3,-2
y = [c_true*np.exp(a_true*i)+d_true*np.exp(b_true*i) for i in wl]
y = np.mat(y)


def Func(abcd:np.ndarray, iput):
    a = abcd[0][0]
    b = abcd[1][0]
    c = abcd[2][0]
    d = abcd[3][0]
    return c*np.exp(a*iput)+d*np.exp(b*iput)


L = Levenberg_Macquardt()
L.Experiment_Wavelength = wl
L.Experiment_Data = y
L.Function = Func
L.Floating_Parameter = np.mat([[-4],[-5],[4],[-4]])
L.Converge()
print(L.Floating_Parameter)