Блог

Вопрос такой:

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

и L, M, N являются целыми константами, 0 <= L, M, N <= 100
, а A, B, C могут быть любыми возможными целыми числами.

Вот приведенные данные:

 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]
 

введите описание изображения здесь

Как использовать градиентный спуск для решения этой тригонометрической функции с несколькими терминами?

Градиентный спуск плохо подходит для оптимизации по целым числам. Вы можете попробовать navie relaxation, где вы решаете с плавающей точкой, и надеюсь, что округленное решение все еще в порядке.

 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()
 

Поскольку проблема будет иметь много локальных минимумов, вам может потребоваться запустить ее несколько раз.

Довольно сложно использовать градиентный спуск для поиска решения этой проблемы, потому что он имеет тенденцию застревать при изменении параметров L, M или N. Градиенты для них могут оттолкнуть его от правильного решения, если только оно уже не очень близко к оптимальному решению.

Есть способы обойти это, такие как basinhopping или случайный поиск, но из-за функции, которую вы пытаетесь изучить, у вас есть лучшая альтернатива.

Поскольку вы пытаетесь изучить функцию синусоиды, вы можете использовать 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)
    )
 

И есть исходная функция!