.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_landweber.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_landweber.py: This is a simulation example of the Landweber class. ==================================================== It is based on the generated supersmooth, smooth and rough signal of Blanchard et al. (2018). .. GENERATED FROM PYTHON SOURCE LINES 7-18 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import EarlyStopping as es from scipy.sparse import dia_matrix import timeit np.random.seed(42) plt.rcParams.update({"font.size": 20}) print('The seed is 42.') .. rst-class:: sphx-glr-script-out .. code-block:: none The seed is 42. .. GENERATED FROM PYTHON SOURCE LINES 19-22 Plot different signals ---------------------- Create diagonal design matrix and supersmooth, smooth and rough signal. Plot the signal. .. GENERATED FROM PYTHON SOURCE LINES 22-42 .. code-block:: Python D = 10000 indices = np.arange(D) + 1 design_matrix = dia_matrix(np.diag(1 / (np.sqrt(indices)))) signal_supersmooth = 5 * np.exp(-0.1 * indices) signal_smooth = 5000 * np.abs(np.sin(0.01 * indices)) * indices ** (-1.6) signal_rough = 250 * np.abs(np.sin(0.002 * indices)) * indices ** (-0.8) plt.figure(figsize=(14, 4)) plt.plot(indices, signal_supersmooth, label="supersmooth signal") plt.plot(indices, signal_smooth, label="smooth signal") plt.plot(indices, signal_rough, label="rough signal") plt.ylabel("Signal") plt.xlabel("Index") plt.xlim([0, 1000]) plt.ylim([0, 1.6]) plt.legend(loc="upper right") plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_001.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 43-46 Generate data and run Landweber ------------------------------- Run the Landweber algorithm and get the early stopping index as well as as the weak/strong balanced oracle. .. GENERATED FROM PYTHON SOURCE LINES 46-103 .. code-block:: Python NOISE_LEVEL = 0.01 noise = np.random.normal(0, NOISE_LEVEL, D) observation_supersmooth = noise + design_matrix @ signal_supersmooth observation_smooth = noise + design_matrix @ signal_smooth observation_rough = noise + design_matrix @ signal_rough models_supersmooth = es.Landweber(design_matrix, observation_supersmooth, true_signal=signal_supersmooth, learning_rate=1, true_noise_level=NOISE_LEVEL) models_smooth = es.Landweber(design_matrix, observation_smooth, true_signal=signal_smooth, learning_rate=1, true_noise_level=NOISE_LEVEL) models_rough = es.Landweber(design_matrix, observation_rough, true_signal=signal_rough, learning_rate=1, true_noise_level=NOISE_LEVEL) # models_supersmooth = es.Landweber( # design_matrix, observation_supersmooth, true_noise_level=NOISE_LEVEL, true_signal=signal_supersmooth # ) # models_smooth = es.Landweber(design_matrix, observation_smooth, true_noise_level=NOISE_LEVEL, true_signal=signal_smooth) # models_rough = es.Landweber(design_matrix, observation_rough, true_noise_level=NOISE_LEVEL, true_signal=signal_rough) max_iteration = 1500 start = timeit.default_timer() models_supersmooth.iterate(max_iteration) stop = timeit.default_timer() print("Time supersmooth: ", stop - start) #models_supersmooth.landweber_estimate_collect start = timeit.default_timer() models_smooth.iterate(max_iteration) stop = timeit.default_timer() print("Time smooth: ", stop - start) start = timeit.default_timer() models_rough.iterate(max_iteration) stop = timeit.default_timer() print("Time rough: ", stop - start) # Stopping index supersmooth_m = models_supersmooth.get_discrepancy_stop(D*(NOISE_LEVEL**2), max_iteration) smooth_m = models_smooth.get_discrepancy_stop(D*(NOISE_LEVEL**2), max_iteration) rough_m = models_rough.get_discrepancy_stop(D*(NOISE_LEVEL**2), max_iteration) # Weak balanced oracle supersmooth_weak_oracle = models_supersmooth.get_weak_balanced_oracle(max_iteration) smooth_weak_oracle = models_smooth.get_weak_balanced_oracle(max_iteration) rough_weak_oracle = models_rough.get_weak_balanced_oracle(max_iteration) # Strong balanced oracle supersmooth_strong_oracle = models_supersmooth.get_strong_balanced_oracle(max_iteration) smooth_strong_oracle = models_smooth.get_strong_balanced_oracle(max_iteration) rough_strong_oracle = models_rough.get_strong_balanced_oracle(max_iteration) .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/EarlyStopping/landweber.py:115: UserWarning: No initial_value is given, using zero by default. warnings.warn("No initial_value is given, using zero by default.", category=UserWarning) /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:603: SparseEfficiencyWarning: splu converted its input to CSC format return splu(A).solve /opt/hostedtoolcache/Python/3.12.7/x64/lib/python3.12/site-packages/scipy/sparse/linalg/_matfuncs.py:76: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format Ainv = spsolve(A, I) Time supersmooth: 3.809562239999991 Time smooth: 3.844058907999994 Time rough: 3.8659052120000013 .. GENERATED FROM PYTHON SOURCE LINES 104-107 Bias/variance decomposition for supersmooth signal -------------------------------------------------- Plot the residuals, weak and strong quantities for the supersmooth signal. .. GENERATED FROM PYTHON SOURCE LINES 107-146 .. code-block:: Python # plt.rcParams.update({'font.size': 18}) # Update the font size fig, axs = plt.subplots(3, 1, figsize=(14, 12)) axs[0].plot(range(0, max_iteration + 1), models_supersmooth.residuals) axs[0].axvline(x=supersmooth_m, color="red", linestyle="--") axs[0].set_xlim([0, 50]) axs[0].set_ylim([0, 20]) axs[0].set_xlabel("Iteration") axs[0].set_ylabel("Residuals") axs[1].plot(range(0, max_iteration + 1), models_supersmooth.strong_risk, color="orange", label="Error") axs[1].plot(range(0, max_iteration + 1), models_supersmooth.strong_bias2, label="$Bias^2$", color="grey") axs[1].plot(range(0, max_iteration + 1), models_supersmooth.strong_variance, label="Variance", color="blue") axs[1].axvline(x=supersmooth_m, color="red", linestyle="--") axs[1].axvline(x=supersmooth_strong_oracle, color="green", linestyle="--") axs[1].set_xlim([0, 50]) axs[1].set_ylim([0, 1]) axs[1].set_xlabel("Iteration") axs[1].set_ylabel("Strong Quantities") axs[2].plot(range(0, max_iteration + 1), models_supersmooth.weak_risk, color="orange", label="Error") axs[2].plot(range(0, max_iteration + 1), models_supersmooth.weak_bias2, label="$Bias^2$", color="grey") axs[2].plot(range(0, max_iteration + 1), models_supersmooth.weak_variance, label="Variance", color="blue") axs[2].axvline(x=supersmooth_m, color="red", linestyle="--", label=r"$\tau$") axs[2].axvline(x=supersmooth_weak_oracle, color="green", linestyle="--", label="$t$ (oracle)") axs[2].set_xlim([0, 400]) axs[2].set_ylim([0, 0.02]) axs[2].set_xlabel("Iteration") axs[2].set_ylabel("Weak Quantities") axs[2].legend() plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_002.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 147-150 Bias/variance decomposition for smooth signal --------------------------------------------- Plot the residuals, weak and strong quantities for the smooth signal. .. GENERATED FROM PYTHON SOURCE LINES 150-186 .. code-block:: Python fig, axs = plt.subplots(3, 1, figsize=(14, 12)) axs[0].plot(range(0, max_iteration + 1), models_smooth.residuals) axs[0].axvline(x=smooth_m, color="red", linestyle="--") axs[0].set_xlim([0, 500]) axs[0].set_ylim([0, 30]) axs[0].set_xlabel("Iteration") axs[0].set_ylabel("Residuals") axs[1].plot(range(0, max_iteration + 1), models_smooth.strong_risk, color="orange", label="Error") axs[1].plot(range(0, max_iteration + 1), models_smooth.strong_bias2, label="$Bias^2$", color="grey") axs[1].plot(range(0, max_iteration + 1), models_smooth.strong_variance, label="Variance", color="blue") axs[1].axvline(x=smooth_m, color="red", linestyle="--") axs[1].axvline(x=smooth_strong_oracle, color="green", linestyle="--") axs[1].set_xlim([0, 500]) axs[1].set_ylim([0, 50]) axs[1].set_xlabel("Iteration") axs[1].set_ylabel("Strong Quantities") axs[2].plot(range(0, max_iteration + 1), models_smooth.weak_risk, color="orange", label="Error") axs[2].plot(range(0, max_iteration + 1), models_smooth.weak_bias2, label="$Bias^2$", color="grey") axs[2].plot(range(0, max_iteration + 1), models_smooth.weak_variance, label="Variance", color="blue") axs[2].axvline(x=smooth_m, color="red", linestyle="--", label=r"$\tau$") axs[2].axvline(x=smooth_weak_oracle, color="green", linestyle="--", label="$t$ (oracle)") axs[2].set_xlim([0, 500]) axs[2].set_ylim([0, 0.5]) axs[2].set_xlabel("Iteration") axs[2].set_ylabel("Weak Quantities") axs[2].legend() plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_003.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 187-190 Bias/variance decomposition for rough signal -------------------------------------------- Plot the residuals, weak and strong quantities for the rough signal. .. GENERATED FROM PYTHON SOURCE LINES 190-227 .. code-block:: Python fig, axs = plt.subplots(3, 1, figsize=(14, 12)) axs[0].plot(range(0, max_iteration + 1), models_rough.residuals) axs[0].axvline(x=rough_m, color="red", linestyle="--") axs[0].set_xlabel("Iteration") axs[0].set_ylabel("Residuals") axs[1].plot(range(0, max_iteration + 1), models_rough.strong_risk, color="orange", label="Error") axs[1].plot(range(0, max_iteration + 1), models_rough.strong_bias2, label="$Bias^2$", color="grey") axs[1].plot(range(0, max_iteration + 1), models_rough.strong_variance, label="Variance", color="blue") axs[1].axvline(x=rough_m, color="red", linestyle="--") axs[1].axvline(x=rough_strong_oracle, color="green", linestyle="--") axs[1].set_ylim([0, 500]) axs[1].set_xlabel("Iteration") axs[1].set_ylabel("Strong Quantities") axs[2].plot(range(0, max_iteration + 1), models_rough.weak_risk, color="orange", label="Error") axs[2].plot(range(0, max_iteration + 1), models_rough.weak_bias2, label="$Bias^2$", color="grey") axs[2].plot(range(0, max_iteration + 1), models_rough.weak_variance, label="Variance", color="blue") axs[2].axvline(x=rough_m, color="red", linestyle="--", label=r"$\tau$") axs[2].axvline(x=rough_weak_oracle, color="green", linestyle="--", label="$t$ (oracle)") axs[2].set_xlim([0, 1200 + 1]) axs[2].set_ylim([0, 1]) axs[2].set_xlabel("Iteration") axs[2].set_ylabel("Weak Quantities") axs[2].legend() plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_004.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 228-232 Bias/variance decomposition for the gravity example --------------------------------------------------- Gravity test problem from the regtools toolbox, see `Hansen (2008) `_ for details. Plot the residuals, weak and strong quantities. .. GENERATED FROM PYTHON SOURCE LINES 232-317 .. code-block:: Python sample_size = 100 # 2**9 a = 0 b = 1 d = 0.25 # Parameter controlling the ill-posedness: the larger, the more ill-posed, default in regtools: d = 0.25 t = (np.arange(1, sample_size + 1) - 0.5) / sample_size s = ((np.arange(1, sample_size + 1) - 0.5) * (b - a)) / sample_size T, S = np.meshgrid(t, s) design = (1 / sample_size) * d * (d**2 * np.ones((sample_size, sample_size)) + (S - T) ** 2) ** (-(3 / 2)) signal = np.sin(np.pi * t) + 0.5 * np.sin(2 * np.pi * t) design_times_signal = design @ signal # Set parameters parameter_size = sample_size max_iteration = 2000 noise_level = 10 ** (-2) # critical_value = sample_size * (noise_level**2) #eigen_values = np.linalg.eig(design) #print(f"The eigenvalues are given by \n {eigen_values}") # Specify number of Monte Carlo runs NUMBER_RUNS = 1 # Create observations noise = np.random.normal(0, noise_level, (sample_size, NUMBER_RUNS)) observation = noise + design_times_signal[:, None] model_gravity = es.Landweber(design, observation[:, 0], learning_rate=1 / 30, true_signal=signal, true_noise_level=noise_level) model_gravity.iterate(max_iteration) # Stopping index m_gravity = model_gravity.get_discrepancy_stop(sample_size * (noise_level**2), max_iteration) print(m_gravity) # Weak balanced oracle weak_oracle_gravity = model_gravity.get_weak_balanced_oracle(max_iteration) # Strong balanced oracle strong_oracle_gravity = model_gravity.get_strong_balanced_oracle(max_iteration) fig, axs = plt.subplots(3, 1, figsize=(14, 12)) axs[0].plot(range(0, max_iteration + 1), model_gravity.residuals) #axs[0].axvline(x=m, color="red", linestyle="--") axs[0].set_xlim([0, 50]) axs[0].set_ylim([0, 20]) axs[0].set_xlabel("Iteration") axs[0].set_ylabel("Residuals") axs[1].plot(range(0, max_iteration + 1), model_gravity.strong_risk, color="orange", label="Error") axs[1].plot(range(0, max_iteration + 1), model_gravity.strong_bias2, label="$Bias^2$", color="grey") axs[1].plot(range(0, max_iteration + 1), model_gravity.strong_variance, label="Variance", color="blue") axs[1].axvline(x=m_gravity, color="red", linestyle="--") axs[1].axvline(x=strong_oracle_gravity, color="green", linestyle="--") #axs[1].set_xlim([0, 50]) axs[1].set_ylim([0, 0.2]) axs[1].set_xlabel("Iteration") axs[1].set_ylabel("Strong Quantities") axs[2].plot(range(0, max_iteration + 1), model_gravity.weak_risk, color="orange", label="Error") axs[2].plot(range(0, max_iteration + 1), model_gravity.weak_bias2, label="$Bias^2$", color="grey") axs[2].plot(range(0, max_iteration + 1), model_gravity.weak_variance, label="Variance", color="blue") axs[2].axvline(x=m_gravity, color="red", linestyle="--", label=r"$\tau$") axs[2].axvline(x=weak_oracle_gravity, color="green", linestyle="--", label="$t$ (oracle)") #axs[2].set_xlim([0, 400]) axs[2].set_ylim([0, 0.002]) axs[2].set_xlabel("Iteration") axs[2].set_ylabel("Weak Quantities") axs[2].legend(loc = "upper right") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_005.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 390 .. GENERATED FROM PYTHON SOURCE LINES 318-321 Bias/variance decomposition for a pertubated diagonal matrix ------------------------------------------------------------ Plot the residuals, weak and strong quantities. .. GENERATED FROM PYTHON SOURCE LINES 321-388 .. code-block:: Python D = 1000 normal_matrix = np.random.normal(0, 0.1, size=(D, D)) indices = np.arange(D) + 1 indices = np.arange(1, D + 1) diagonal_values = 1 / np.sqrt(indices) np.fill_diagonal(normal_matrix, diagonal_values) design_matrix = normal_matrix NOISE_LEVEL = 0.1 noise = np.random.normal(0, NOISE_LEVEL, D) indices = np.arange(D) + 1 signal_supersmooth = 5 * np.exp(-0.1 * indices) response = noise + design_matrix @ signal_supersmooth max_iteration = 250 model_pertubation = es.Landweber(design_matrix, response, learning_rate = 0.01 , true_signal=signal_supersmooth, true_noise_level=NOISE_LEVEL) model_pertubation.iterate(max_iteration) # Stopping index m_pertubation = model_pertubation.get_discrepancy_stop(D*(NOISE_LEVEL**2), max_iteration) # Weak balanced oracle weak_oracle_pertubation = model_pertubation.get_weak_balanced_oracle(max_iteration) # Strong balanced oracle strong_oracle_pertubation = model_pertubation.get_strong_balanced_oracle(max_iteration) fig, axs = plt.subplots(3, 1, figsize=(14, 12)) axs[0].plot(range(0, max_iteration + 1), model_pertubation.residuals) # axs[0].axvline(x=m, color="red", linestyle="--") # axs[0].set_xlim([0, 50]) axs[0].set_ylim([0, 20]) axs[0].set_xlabel("Iteration") axs[0].set_ylabel("Residuals") axs[1].plot(range(0, max_iteration + 1), model_pertubation.strong_risk, color="orange", label="Error") axs[1].plot(range(0, max_iteration + 1), model_pertubation.strong_bias2, label="$Bias^2$", color="grey") axs[1].plot(range(0, max_iteration + 1), model_pertubation.strong_variance, label="Variance", color="blue") axs[1].axvline(x=m_pertubation, color="red", linestyle="--") axs[1].axvline(x=strong_oracle_pertubation, color="green", linestyle="--") # axs[1].set_xlim([0, 50]) axs[1].set_ylim([0, 400]) axs[1].set_xlabel("Iteration") axs[1].set_ylabel("Strong Quantities") axs[2].plot(range(0, max_iteration + 1), model_pertubation.weak_risk, color="orange", label="Error") axs[2].plot(range(0, max_iteration + 1), model_pertubation.weak_bias2, label="$Bias^2$", color="grey") axs[2].plot(range(0, max_iteration + 1), model_pertubation.weak_variance, label="Variance", color="blue") axs[2].axvline(x=m_pertubation, color="red", linestyle="--", label=r"$\tau$") axs[2].axvline(x=weak_oracle_pertubation, color="green", linestyle="--", label="$t$ (oracle)") # axs[2].set_xlim([0, 400]) axs[2].set_ylim([0, 50]) axs[2].set_xlabel("Iteration") axs[2].set_ylabel("Weak Quantities") axs[2].legend(loc = "upper right") plt.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_landweber_006.png :alt: plot landweber :srcset: /auto_examples/images/sphx_glr_plot_landweber_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 54.372 seconds) .. _sphx_glr_download_auto_examples_plot_landweber.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_landweber.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_landweber.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_landweber.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_