From 3dd1acb6432fbd903f08ba288a43dacdff723209 Mon Sep 17 00:00:00 2001 From: Jan Wilhelm Date: Tue, 30 May 2023 22:45:32 +0200 Subject: [PATCH 1/2] Introduce regularization --- GX-TimeFrequency/src/minimax_grids.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GX-TimeFrequency/src/minimax_grids.F90 b/GX-TimeFrequency/src/minimax_grids.F90 index e94de9aa..7f9c1ed8 100644 --- a/GX-TimeFrequency/src/minimax_grids.F90 +++ b/GX-TimeFrequency/src/minimax_grids.F90 @@ -253,7 +253,7 @@ subroutine get_transformation_weights(num_points, tau_points, omega_points, weig integer :: i_node, i_point, j_point, k_point, & num_x_nodes integer, parameter :: nodes_factor = 200 - real(kind=dp) :: current_point, x_factor + real(kind=dp) :: current_point, x_factor, regularization real(kind=dp), allocatable, dimension(:) :: weights_work, x_mu, psi real(kind=dp), allocatable, dimension(:, :) :: mat_A @@ -316,9 +316,10 @@ subroutine get_transformation_weights(num_points, tau_points, omega_points, weig ! integration weights = (V Sigma^-1 U^T)*psi ! 1) V * Sigma^-1 + regularization = 0.01_dp do j_point = 1, num_points do k_point = 1, num_points - mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)/vec_S(j_point) + mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)/(vec_S(j_point)+regularization) end do ! k_point end do ! j_point From bbbef1ba0e26c06370671693ae445231780a31f0 Mon Sep 17 00:00:00 2001 From: Jan Wilhelm Date: Tue, 30 May 2023 22:51:51 +0200 Subject: [PATCH 2/2] correct for Wiener filter https://en.wikipedia.org/wiki/Ridge_regression See "Relation to singular-value decomposition and Wiener filter" --- GX-TimeFrequency/src/minimax_grids.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/GX-TimeFrequency/src/minimax_grids.F90 b/GX-TimeFrequency/src/minimax_grids.F90 index 7f9c1ed8..3ff09710 100644 --- a/GX-TimeFrequency/src/minimax_grids.F90 +++ b/GX-TimeFrequency/src/minimax_grids.F90 @@ -319,7 +319,8 @@ subroutine get_transformation_weights(num_points, tau_points, omega_points, weig regularization = 0.01_dp do j_point = 1, num_points do k_point = 1, num_points - mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)/(vec_S(j_point)+regularization) + mat_VT_s(k_point, j_point) = mat_VT(j_point, k_point)*vec_S(j_point) / & + (vec_S(j_point)**2+regularization**2) end do ! k_point end do ! j_point