如何使用梯度下降来解决这个多项三角函数？

Question

问题是这样的：

f(x) = A sin(2π * L * x) + B cos(2π * M * x) + C sin(2π * N * x)

和 L,M,N 是常量整数，0

这是给定的数据：

x = [0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99]
y = [4,1.240062433,-0.7829654986,-1.332487982,-0.3337640721,1.618033989,3.512512389,4.341307895,3.515268061,1.118929599,-2.097886967,-4.990538967,-6.450324073,-5.831575611,-3.211486891,0.6180339887,4.425660706,6.980842552,7.493970785,5.891593744,2.824429495,-0.5926374511,-3.207870455,-4.263694544,-3.667432785,-2,-0.2617162175,0.5445886005,-0.169441247,-2.323237059,-5.175570505,-7.59471091,-8.488730333,-7.23200463,-3.924327772,0.6180339887,5.138501587,8.38127157,9.532377045,8.495765687,5.902113033,2.849529206,0.4768388529,-0.46697525,0.106795821,1.618033989,3.071952496,3.475795162,2.255463709,-0.4905371745,-4,-7.117914956,-8.727599664,-8.178077181,-5.544088451,-1.618033989,2.365340134,5.169257268,5.995297102,4.758922924,2.097886967,-0.8873135564,-3.06024109,-3.678989552,-2.666365632,-0.6180339887,1.452191817,2.529722611,2.016594378,-0.01374122059,-2.824429495,-5.285215072,-6.302694708,-5.246870619,-2.210419738,2,6.13956874,8.965976562,9.68000641,8.201089581,5.175570505,1.716858387,-1.02183483,-2.278560533,-1.953524751,-0.6180339887,0.7393509358,1.129293593,-0.02181188158,-2.617913164,-5.902113033,-8.727381729,-9.987404016,-9.043589913,-5.984648344,-1.618033989,2.805900027,6.034770001,7.255101454,6.368389697]

enter image description here

如何使用梯度下降来求解这个多项三角函数？

Answer 1

梯度下降不太适合对整数进行优化。您可以尝试在浮点数中求解的 navie 放松，并希望圆形求解仍然可以。

from autograd import grad, numpy as jnp
import numpy as np

def cast(params):
    [A, B, C, L, M, N] = params
    L = jnp.minimum(jnp.abs(L), 100)
    M = jnp.minimum(jnp.abs(M), 100)
    N = jnp.minimum(jnp.abs(N), 100)
    return A, B, C, L, M, N

def pred(params, x):
  [A, B, C, L, M, N] = cast(params)
  return A *jnp.sin(2 * jnp.pi * L * x) + B*jnp.cos(2*jnp.pi * M * x) + C * jnp.sin(2 * jnp.pi * N * x)

x = [0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99]
y = [4,1.240062433,-0.7829654986,-1.332487982,-0.3337640721,1.618033989,3.512512389,4.341307895,3.515268061,1.118929599,-2.097886967,-4.990538967,-6.450324073,-5.831575611,-3.211486891,0.6180339887,4.425660706,6.980842552,7.493970785,5.891593744,2.824429495,-0.5926374511,-3.207870455,-4.263694544,-3.667432785,-2,-0.2617162175,0.5445886005,-0.169441247,-2.323237059,-5.175570505,-7.59471091,-8.488730333,-7.23200463,-3.924327772,0.6180339887,5.138501587,8.38127157,9.532377045,8.495765687,5.902113033,2.849529206,0.4768388529,-0.46697525,0.106795821,1.618033989,3.071952496,3.475795162,2.255463709,-0.4905371745,-4,-7.117914956,-8.727599664,-8.178077181,-5.544088451,-1.618033989,2.365340134,5.169257268,5.995297102,4.758922924,2.097886967,-0.8873135564,-3.06024109,-3.678989552,-2.666365632,-0.6180339887,1.452191817,2.529722611,2.016594378,-0.01374122059,-2.824429495,-5.285215072,-6.302694708,-5.246870619,-2.210419738,2,6.13956874,8.965976562,9.68000641,8.201089581,5.175570505,1.716858387,-1.02183483,-2.278560533,-1.953524751,-0.6180339887,0.7393509358,1.129293593,-0.02181188158,-2.617913164,-5.902113033,-8.727381729,-9.987404016,-9.043589913,-5.984648344,-1.618033989,2.805900027,6.034770001,7.255101454,6.368389697]

def loss(params):
  p = pred(params, np.array(x))
  return jnp.mean((np.array(y)-p)**2)

params = np.array([np.random.random()*100 for _ in range(6)])
for _ in range(10000):
  g = grad(loss)

  params = params - 0.001*g(params)
  print("Relaxed solution", cast(params), "loss=", loss(params))

  constrained_params = np.round(cast(params))
  print("Integer solution", constrained_params, "loss=", loss(constrained_params))
  print()

由于问题会有很多局部最小值，您可能需要多次运行它。

Answer 2

使用梯度下降来解决这个问题是相当困难的，因为它在改变 L、M 或 N 参数时往往会卡住。那些梯度可能会将其推离正确的解决方案，除非它已经非常接近最优解决方案。

有一些方法可以解决这个问题，例如盆地跳跃或随机搜索，但由于您尝试学习的功能，您有更好的选择。

由于您正在尝试学习正弦函数，您可以使用FFT 来查找正弦波的频率。一旦有了这些频率，您就可以找到用于生成相同正弦波的幅度和相位。

请原谅这段代码的混乱，这是我第一次使用 FFT。

import scipy.fft
import numpy as np
import math
import matplotlib.pyplot as plt

def get_top_frequencies(x, y, num_freqs):
    x = np.array(x)
    y = np.array(y)

    # Find timestep (assume constant timestep)
    dt = abs(x[0] - x[-1]) / (len(x) - 1)
    
    # Take discrete FFT of y
    spectral = scipy.fft.fft(y)
    freq = scipy.fft.fftfreq(y.shape[0], d=dt)

    # Cut off top half of frequencies. Assumes input signal is real, and not complex.
    spectral = spectral[:int(spectral.shape[0] / 2)]
    # Double amplitudes to correct for cutting off top half.
    spectral *= 2
    
    # Adjust amplitude by sampling timestep
    spectral *= dt
    
    # Get ampitudes for all frequencies. This is taking the magnitude of the complex number
    spectral_amplitude = np.abs(spectral)
    
    # Pick frequencies with highest amplitudes
    highest_idx = np.argsort(spectral_amplitude)[::-1][:num_freqs]
    
    # Find amplitude, frequency, and phase components of each term
    highest_amplitude = spectral_amplitude[highest_idx]
    highest_freq = freq[highest_idx]
    highest_phase = np.angle(spectral[highest_idx]) / math.pi

    # Convert it into a Python function
    function = ["def func(x):", "return ("]
    
    for i, components in enumerate(zip(highest_amplitude, highest_freq, highest_phase)):
        amplitude, freq, phase = components
        plus_sign = " +" if i != (num_freqs - 1) else ""
        term = f"{amplitude:.2f} * math.cos(2 * math.pi * {freq:.2f} * x + math.pi * {phase:.2f}){plus_sign}"
        function.append("    " + term)
    function.append(")")
    return "\n    ".join(function)


x = [0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99]
y = [4,1.240062433,-0.7829654986,-1.332487982,-0.3337640721,1.618033989,3.512512389,4.341307895,3.515268061,1.118929599,-2.097886967,-4.990538967,-6.450324073,-5.831575611,-3.211486891,0.6180339887,4.425660706,6.980842552,7.493970785,5.891593744,2.824429495,-0.5926374511,-3.207870455,-4.263694544,-3.667432785,-2,-0.2617162175,0.5445886005,-0.169441247,-2.323237059,-5.175570505,-7.59471091,-8.488730333,-7.23200463,-3.924327772,0.6180339887,5.138501587,8.38127157,9.532377045,8.495765687,5.902113033,2.849529206,0.4768388529,-0.46697525,0.106795821,1.618033989,3.071952496,3.475795162,2.255463709,-0.4905371745,-4,-7.117914956,-8.727599664,-8.178077181,-5.544088451,-1.618033989,2.365340134,5.169257268,5.995297102,4.758922924,2.097886967,-0.8873135564,-3.06024109,-3.678989552,-2.666365632,-0.6180339887,1.452191817,2.529722611,2.016594378,-0.01374122059,-2.824429495,-5.285215072,-6.302694708,-5.246870619,-2.210419738,2,6.13956874,8.965976562,9.68000641,8.201089581,5.175570505,1.716858387,-1.02183483,-2.278560533,-1.953524751,-0.6180339887,0.7393509358,1.129293593,-0.02181188158,-2.617913164,-5.902113033,-8.727381729,-9.987404016,-9.043589913,-5.984648344,-1.618033989,2.805900027,6.034770001,7.255101454,6.368389697]
print(get_top_frequencies(x, y, 3))

产生这个函数：

def func(x):
    return (
        5.00 * math.cos(2 * math.pi * 10.00 * x + math.pi * 0.50) +
        4.00 * math.cos(2 * math.pi * 5.00 * x + math.pi * -0.00) +
        2.00 * math.cos(2 * math.pi * 3.00 * x + math.pi * -0.50)
    )

这不是您指定的格式 - 您要求两个 sin 和一个 cos，并且没有相位参数。但是，使用三角恒等式 cos(x) = sin(pi/2 - x)，您可以将其转换为与您想要的匹配的等效表达式：

def func(x):
    return (
        5.00 * math.sin(2 * math.pi * -10.00 * x) +
        4.00 * math.cos(2 * math.pi * 5.00 * x) +
        2.00 * math.sin(2 * math.pi * 3.00 * x)
    )

还有原来的功能！