From ef16dc848dd808ffb475fdd0ee5281aae930ca78 Mon Sep 17 00:00:00 2001 From: IgorKoh Date: Fri, 9 Apr 2021 15:57:00 +0200 Subject: [PATCH] Partial implementation of stevengj suggestions --- Project.toml | 11 ++ src/Interpolate.jl | 351 ++---------------------------------- src/NormalHermiteSplines.jl | 3 +- src/ReproducingKernels.jl | 8 +- src/Utils.jl | 331 ++++++++++++++++++++++++++++++++++ 5 files changed, 359 insertions(+), 345 deletions(-) create mode 100644 src/Utils.jl diff --git a/Project.toml b/Project.toml index 46cc7bb1..b77b5dfe 100644 --- a/Project.toml +++ b/Project.toml @@ -5,7 +5,18 @@ authors = ["Igor Kohanovsky "] version = "0.5.1" [deps] +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" +Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004" +IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] diff --git a/src/Interpolate.jl b/src/Interpolate.jl index a4ca9ad8..097300ca 100644 --- a/src/Interpolate.jl +++ b/src/Interpolate.jl @@ -1,288 +1,3 @@ -function _estimate_accuracy(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0} - n = size(spline._nodes, 1) - m = size(spline._nodes, 2) - nodes = similar(spline._nodes) - @inbounds for i = 1:n - for j = 1:m - nodes[i,j] = spline._min_bound[i] + spline._compression * spline._nodes[i,j] - end - end - σ = evaluate(spline, nodes) - # calculating a value of the Relative Maximum Absolute Error (RMAE) of interpolation - # at the function value interpolation nodes. - fun_max = maximum(abs.(spline._values)) - if fun_max > 0.0 - rmae = maximum(abs.(spline._values .- σ)) / fun_max - else - rmae = maximum(abs.(spline._values .- σ)) - end - rmae = rmae > eps(T) ? rmae : eps(T) - res = -floor(log10(rmae)) - 1 - if res <= 0 - res = 0 - end - return trunc(Int, res) -end - -function _estimate_ε(nodes::Matrix{T}) where T <: AbstractFloat - n = size(nodes, 1) - n_1 = size(nodes, 2) - ε = T(0.0) - @inbounds for i = 1:n_1 - for j = i:n_1 - ε += norm(nodes[:,i] .- nodes[:,j]) - end - end - if ε > T(0.0) - ε *= T(n)^(1.0/n) / T(n_1)^(5.0/3.0) - else - ε = T(1.0) - end - return ε -end - -function _estimate_ε(nodes::Matrix{T}, - d_nodes::Matrix{T} - ) where T <: AbstractFloat - return _estimate_ε([nodes 0.1 .* d_nodes]) -end - -function _estimate_epsilon(nodes::Matrix{T}, - kernel::RK) where {T <: AbstractFloat, RK <: ReproducingKernel_0} - n = size(nodes, 1) - n_1 = size(nodes, 2) - min_bound = Vector{T}(undef, n) - compression::T = 0 - @inbounds for i = 1:n - min_bound[i] = nodes[i,1] - maxx::T = nodes[i,1] - for j = 2:n_1 - min_bound[i] = min(min_bound[i], nodes[i,j]) - maxx = max(maxx, nodes[i,j]) - end - compression = max(compression, maxx - min_bound[i]) - end - - if compression <= eps(T(1.0)) - error("Cannot estimate_epsilon: `nodes` data are not correct.") - end - - t_nodes = similar(nodes) - @inbounds for i = 1:n - for j = 1:n_1 - t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression - end - end - ε = _estimate_ε(t_nodes) - if isa(kernel, RK_H1) - ε *= T(1.5) - elseif isa(kernel, RK_H2) - ε *= T(2.0) - end - return ε -end - -function _estimate_epsilon(nodes::Matrix{T}, - d_nodes::Matrix{T}, - kernel::RK) where {T <: AbstractFloat, RK <: ReproducingKernel_1} - n = size(nodes, 1) - n_1 = size(nodes, 2) - n_2 = size(d_nodes, 2) - - min_bound = Vector{T}(undef, n) - compression::T = 0 - @inbounds for i = 1:n - min_bound[i] = nodes[i,1] - maxx::T = nodes[i,1] - for j = 2:n_1 - min_bound[i] = min(min_bound[i], nodes[i,j]) - maxx = max(maxx, nodes[i,j]) - end - for j = 1:n_2 - min_bound[i] = min(min_bound[i], d_nodes[i,j]) - maxx = max(maxx, d_nodes[i,j]) - end - compression = max(compression, maxx - min_bound[i]) - end - - if compression <= eps(T(1.0)) - error("Cannot estimate_epsilon: `nodes`, `d_nodes` data are not correct.") - end - - t_nodes = similar(nodes) - t_d_nodes = similar(d_nodes) - @inbounds for i = 1:n - for j = 1:n_1 - t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression - end - for j = 1:n_2 - t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression - end - end - - ε = _estimate_ε(t_nodes, t_d_nodes) - if isa(kernel, RK_H1) - ε *= T(2.0) - elseif isa(kernel, RK_H2) - ε *= T(2.5) - end - - return ε -end - -function _get_gram(nodes::Matrix{T}, - kernel::RK - ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} - n = size(nodes, 1) - n_1 = size(nodes, 2) - min_bound = Vector{T}(undef, n) - compression::T = 0 - @inbounds for i = 1:n - min_bound[i] = nodes[i,1] - maxx::T = nodes[i,1] - for j = 2:n_1 - min_bound[i] = min(min_bound[i], nodes[i,j]) - maxx = max(maxx, nodes[i,j]) - end - compression = max(compression, maxx - min_bound[i]) - end - - if compression <= eps(T(1.0)) - error("Cannot prepare the spline: `nodes` data are not correct.") - end - - t_nodes = similar(nodes) - @inbounds for i = 1:n - for j = 1:n_1 - t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression - end - end - - if T(kernel.ε) == T(0.0) - ε = _estimate_ε(t_nodes) - if isa(kernel, RK_H0) - kernel = RK_H0(ε) - elseif isa(kernel, RK_H1) - ε *= T(1.5) - kernel = RK_H1(ε) - elseif isa(kernel, RK_H2) - ε *= T(2.0) - kernel = RK_H2(ε) - else - error("incorrect `kernel` type.") - end - end - - return _gram(t_nodes, kernel) -end - -function _get_gram(nodes::Matrix{T}, - d_nodes::Matrix{T}, - es::Matrix{T}, - kernel::RK - ) where {T <: AbstractFloat, RK <: ReproducingKernel_1} - n = size(nodes, 1) - n_1 = size(nodes, 2) - n_2 = size(d_nodes, 2) - - if(size(es, 2) != n_2) - error("Number of derivative directions does not correspond to the number of derivative nodes.") - end - - t_es = similar(es) - try - @inbounds for i = 1:n_2 - t_es[:,i] = es[:,i] ./ norm(es[:,i]) - end - catch - error("Cannot normalize derivative direction: zero direction vector.") - end - - min_bound = Vector{T}(undef, n) - compression::T = 0 - @inbounds for i = 1:n - min_bound[i] = nodes[i,1] - maxx::T = nodes[i,1] - for j = 2:n_1 - min_bound[i] = min(min_bound[i], nodes[i,j]) - maxx = max(maxx, nodes[i,j]) - end - for j = 1:n_2 - min_bound[i] = min(min_bound[i], d_nodes[i,j]) - maxx = max(maxx, d_nodes[i,j]) - end - compression = max(compression, maxx - min_bound[i]) - end - - if compression <= eps(T(1.0)) - error("Cannot prepare the spline: `nodes`, `d_nodes` data are not correct.") - end - - t_nodes = similar(nodes) - t_d_nodes = similar(d_nodes) - @inbounds for i = 1:n - for j = 1:n_1 - t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression - end - for j = 1:n_2 - t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression - end - end - - if T(kernel.ε) == T(0.0) - ε = _estimate_ε(t_nodes, t_d_nodes) - if isa(kernel, RK_H1) - ε *= T(2.0) - kernel = RK_H1(ε) - elseif isa(kernel, RK_H2) - ε *= T(2.5) - kernel = RK_H2(ε) - else - error("incorrect `kernel` type.") - end - end - - return _gram(t_nodes, t_d_nodes, t_es, kernel) -end - -function _get_cond(nodes::Matrix{T}, - kernel::RK - ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} - cond = T(0.0) - gram = _get_gram(nodes, kernel) - try - evs = svdvals!(gram) - maxevs = maximum(evs) - minevs = minimum(evs) - if minevs > T(0.0) - cond = maxevs / minevs - cond = 10.0^floor(log10(cond)) - end - catch - end - return cond -end - -function _get_cond(nodes::Matrix{T}, - d_nodes::Matrix{T}, - es::Matrix{T}, - kernel::RK - ) where {T <: AbstractFloat, RK <: ReproducingKernel_1} - cond = T(0.0) - gram = _get_gram(nodes, kernel) - try - evs = svdvals!(gram) - maxevs = maximum(evs) - minevs = minimum(evs) - if minevs > T(0.0) - cond = maxevs / minevs - cond = 10.0^floor(log10(cond)) - end - catch - end - return cond -end - function _prepare(nodes::Matrix{T}, kernel::RK ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} @@ -299,14 +14,13 @@ function _prepare(nodes::Matrix{T}, end compression = max(compression, maxx - min_bound[i]) end - if compression <= eps(T(1.0)) error("Cannot prepare the spline: `nodes` data are not correct.") end t_nodes = similar(nodes) - @inbounds for i = 1:n - for j = 1:n_1 + @inbounds for j = 1:n_1 + for i = 1:n t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression end end @@ -428,11 +142,15 @@ function _prepare(nodes::Matrix{T}, t_nodes = similar(nodes) t_d_nodes = similar(d_nodes) - @inbounds for i = 1:n - for j = 1:n_1 + + @inbounds for j = 1:n_1 + for i = 1:n t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression end - for j = 1:n_2 + end + + @inbounds for j = 1:n_2 + for i = 1:n t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression end end @@ -534,8 +252,8 @@ function _evaluate(spline::NormalSpline{T, RK}, m = size(points, 2) pts = similar(points) - @inbounds for i = 1:n - for j = 1:m + @inbounds for j = 1:m + for i = 1:n pts[i,j] = (points[i,j] - spline._min_bound[i]) / spline._compression end end @@ -671,50 +389,3 @@ function _evaluate_gradient(spline::NormalSpline{T, RK}, end return grad end - -# ``` -# Get estimation of the Gram matrix condition number -# Brás, C.P., Hager, W.W. & Júdice, J.J. An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP 20, 791–809 (2012). -# https://link.springer.com/article/10.1007/s11750-010-0161-9 -# ``` -function _estimate_cond(gram::Matrix{T}, - chol::LinearAlgebra.Cholesky{T,Array{T,2}}, - nit = 3 - ) where T <: AbstractFloat - if isnothing(gram) - throw(DomainError(gram, "Parameter `gram` is `nothing'.")) - end - if isnothing(chol) - throw(DomainError(chol, "Parameter `chol` is `nothing'.")) - end - mat_norm = norm(gram, 1) - n = size(gram, 1) - x = Vector{T}(undef, n) - @. x = T(1.0) / T(n) - z = Vector{T}(undef, n) - gamma = T(0.0) - for it = 1:nit - z = ldiv!(z, chol, x) - gamma = T(0.0); - for i = 1:n - gamma += abs(z[i]) - z[i] = sign(z[i]) - end - z = ldiv!(z, chol, copy(z)) - zx = z ⋅ x - idx = 1 - for i = 1:n - z[i] = abs(z[i]) - if z[i] > z[idx] - idx = i - end - end - if z[idx] <= zx - break - end - @. x = T(0.0) - x[idx] = T(1.0) - end - cond = T(10.0)^floor(log10(mat_norm * gamma)) - return cond -end diff --git a/src/NormalHermiteSplines.jl b/src/NormalHermiteSplines.jl index 802e29a3..99cf4300 100644 --- a/src/NormalHermiteSplines.jl +++ b/src/NormalHermiteSplines.jl @@ -10,7 +10,7 @@ export estimate_accuracy export evaluate_derivative # -- #### -#include("./examples/Main.jl") +include("./examples/Main.jl") using LinearAlgebra @@ -56,6 +56,7 @@ end include("./ReproducingKernels.jl") include("./GramMatrix.jl") +include("./Utils.jl") include("./Interpolate.jl") """ diff --git a/src/ReproducingKernels.jl b/src/ReproducingKernels.jl index 795dc638..f3470ffd 100644 --- a/src/ReproducingKernels.jl +++ b/src/ReproducingKernels.jl @@ -89,7 +89,7 @@ end ξ::Vector{T} ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} defined::Bool = false - x::T = kernel.ε * norm(ξ .- η) + x::T = kernel.ε * sqrt(sum((ξ - η) .^ 2)) if isa(kernel, RK_H2) defined = true value = (T(3.0) + x * (T(3.0) + x)) * exp(-x) @@ -140,7 +140,7 @@ end # Note: Derivative of spline built with reproducing kernel RK_H0 does not exist at the spline nodes. value::T = T(0.0) defined::Bool = false - normt = norm(η .- ξ) + normt = sqrt(sum((η - ξ) .^ 2)) x = kernel.ε * normt if isa(kernel, RK_H2) defined = true @@ -173,7 +173,7 @@ end defined::Bool = false if isa(kernel, RK_H2) defined = true - x = kernel.ε * norm(η .- ξ) + x = kernel.ε * sqrt(sum((η - ξ) .^ 2)) if r == k if x > T(0.0) value = kernel.ε^2 * exp(-x) * (T(1.0) + x - (kernel.ε * (ξ[r] - η[r]))^2) @@ -190,7 +190,7 @@ end end if isa(kernel, RK_H1) defined = true - t = norm(η .- ξ) + t = sqrt(sum((η - ξ) .^ 2)) x = kernel.ε * t if r == k if t > T(0.0) diff --git a/src/Utils.jl b/src/Utils.jl new file mode 100644 index 00000000..135180c8 --- /dev/null +++ b/src/Utils.jl @@ -0,0 +1,331 @@ +function _estimate_accuracy(spline::NormalSpline{T, RK}) where {T <: AbstractFloat, RK <: ReproducingKernel_0} + n = size(spline._nodes, 1) + m = size(spline._nodes, 2) + nodes = similar(spline._nodes) + @inbounds for i = 1:n + for j = 1:m + nodes[i,j] = spline._min_bound[i] + spline._compression * spline._nodes[i,j] + end + end + σ = evaluate(spline, nodes) + # calculating a value of the Relative Maximum Absolute Error (RMAE) of interpolation + # at the function value interpolation nodes. + fun_max = maximum(abs.(spline._values)) + if fun_max > 0.0 + rmae = maximum(abs.(spline._values .- σ)) / fun_max + else + rmae = maximum(abs.(spline._values .- σ)) + end + rmae = rmae > eps(T) ? rmae : eps(T) + res = -floor(log10(rmae)) - 1 + if res <= 0 + res = 0 + end + return trunc(Int, res) +end + +function _estimate_ε(nodes::Matrix{T}) where T <: AbstractFloat + n = size(nodes, 1) + n_1 = size(nodes, 2) + ε = T(0.0) + @inbounds for i = 1:n_1 + for j = i:n_1 + ε += norm(nodes[:,i] .- nodes[:,j]) + end + end + if ε > T(0.0) + ε *= T(n)^(1.0/n) / T(n_1)^(5.0/3.0) + else + ε = T(1.0) + end + return ε +end + +function _estimate_ε(nodes::Matrix{T}, + d_nodes::Matrix{T} + ) where T <: AbstractFloat + return _estimate_ε([nodes 0.1 .* d_nodes]) +end + +function _estimate_epsilon(nodes::Matrix{T}, + kernel::RK) where {T <: AbstractFloat, RK <: ReproducingKernel_0} + n = size(nodes, 1) + n_1 = size(nodes, 2) + min_bound = Vector{T}(undef, n) + compression::T = 0 + @inbounds for i = 1:n + min_bound[i] = nodes[i,1] + maxx::T = nodes[i,1] + for j = 2:n_1 + min_bound[i] = min(min_bound[i], nodes[i,j]) + maxx = max(maxx, nodes[i,j]) + end + compression = max(compression, maxx - min_bound[i]) + end + + if compression <= eps(T(1.0)) + error("Cannot estimate_epsilon: `nodes` data are not correct.") + end + + t_nodes = similar(nodes) + @inbounds for i = 1:n + for j = 1:n_1 + t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression + end + end + ε = _estimate_ε(t_nodes) + if isa(kernel, RK_H1) + ε *= T(1.5) + elseif isa(kernel, RK_H2) + ε *= T(2.0) + end + return ε +end + +function _estimate_epsilon(nodes::Matrix{T}, + d_nodes::Matrix{T}, + kernel::RK) where {T <: AbstractFloat, RK <: ReproducingKernel_1} + n = size(nodes, 1) + n_1 = size(nodes, 2) + n_2 = size(d_nodes, 2) + + min_bound = Vector{T}(undef, n) + compression::T = 0 + @inbounds for i = 1:n + min_bound[i] = nodes[i,1] + maxx::T = nodes[i,1] + for j = 2:n_1 + min_bound[i] = min(min_bound[i], nodes[i,j]) + maxx = max(maxx, nodes[i,j]) + end + for j = 1:n_2 + min_bound[i] = min(min_bound[i], d_nodes[i,j]) + maxx = max(maxx, d_nodes[i,j]) + end + compression = max(compression, maxx - min_bound[i]) + end + + if compression <= eps(T(1.0)) + error("Cannot estimate_epsilon: `nodes`, `d_nodes` data are not correct.") + end + + t_nodes = similar(nodes) + t_d_nodes = similar(d_nodes) + @inbounds for i = 1:n + for j = 1:n_1 + t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression + end + for j = 1:n_2 + t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression + end + end + + ε = _estimate_ε(t_nodes, t_d_nodes) + if isa(kernel, RK_H1) + ε *= T(2.0) + elseif isa(kernel, RK_H2) + ε *= T(2.5) + end + + return ε +end + +function _get_gram(nodes::Matrix{T}, + kernel::RK + ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} + n = size(nodes, 1) + n_1 = size(nodes, 2) + min_bound = Vector{T}(undef, n) + compression::T = 0 + @inbounds for i = 1:n + min_bound[i] = nodes[i,1] + maxx::T = nodes[i,1] + for j = 2:n_1 + min_bound[i] = min(min_bound[i], nodes[i,j]) + maxx = max(maxx, nodes[i,j]) + end + compression = max(compression, maxx - min_bound[i]) + end + + if compression <= eps(T(1.0)) + error("Cannot prepare the spline: `nodes` data are not correct.") + end + + t_nodes = similar(nodes) + @inbounds for i = 1:n + for j = 1:n_1 + t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression + end + end + + if T(kernel.ε) == T(0.0) + ε = _estimate_ε(t_nodes) + if isa(kernel, RK_H0) + kernel = RK_H0(ε) + elseif isa(kernel, RK_H1) + ε *= T(1.5) + kernel = RK_H1(ε) + elseif isa(kernel, RK_H2) + ε *= T(2.0) + kernel = RK_H2(ε) + else + error("incorrect `kernel` type.") + end + end + + return _gram(t_nodes, kernel) +end + +function _get_gram(nodes::Matrix{T}, + d_nodes::Matrix{T}, + es::Matrix{T}, + kernel::RK + ) where {T <: AbstractFloat, RK <: ReproducingKernel_1} + n = size(nodes, 1) + n_1 = size(nodes, 2) + n_2 = size(d_nodes, 2) + + if(size(es, 2) != n_2) + error("Number of derivative directions does not correspond to the number of derivative nodes.") + end + + t_es = similar(es) + try + @inbounds for i = 1:n_2 + t_es[:,i] = es[:,i] ./ norm(es[:,i]) + end + catch + error("Cannot normalize derivative direction: zero direction vector.") + end + + min_bound = Vector{T}(undef, n) + compression::T = 0 + @inbounds for i = 1:n + min_bound[i] = nodes[i,1] + maxx::T = nodes[i,1] + for j = 2:n_1 + min_bound[i] = min(min_bound[i], nodes[i,j]) + maxx = max(maxx, nodes[i,j]) + end + for j = 1:n_2 + min_bound[i] = min(min_bound[i], d_nodes[i,j]) + maxx = max(maxx, d_nodes[i,j]) + end + compression = max(compression, maxx - min_bound[i]) + end + + if compression <= eps(T(1.0)) + error("Cannot prepare the spline: `nodes`, `d_nodes` data are not correct.") + end + + t_nodes = similar(nodes) + t_d_nodes = similar(d_nodes) + @inbounds for i = 1:n + for j = 1:n_1 + t_nodes[i,j] = (nodes[i,j] - min_bound[i]) / compression + end + for j = 1:n_2 + t_d_nodes[i,j] = (d_nodes[i,j] - min_bound[i]) / compression + end + end + + if T(kernel.ε) == T(0.0) + ε = _estimate_ε(t_nodes, t_d_nodes) + if isa(kernel, RK_H1) + ε *= T(2.0) + kernel = RK_H1(ε) + elseif isa(kernel, RK_H2) + ε *= T(2.5) + kernel = RK_H2(ε) + else + error("incorrect `kernel` type.") + end + end + + return _gram(t_nodes, t_d_nodes, t_es, kernel) +end + +function _get_cond(nodes::Matrix{T}, + kernel::RK + ) where {T <: AbstractFloat, RK <: ReproducingKernel_0} + cond = T(0.0) + gram = _get_gram(nodes, kernel) + try + evs = svdvals!(gram) + maxevs = maximum(evs) + minevs = minimum(evs) + if minevs > T(0.0) + cond = maxevs / minevs + cond = 10.0^floor(log10(cond)) + end + catch + end + return cond +end + +function _get_cond(nodes::Matrix{T}, + d_nodes::Matrix{T}, + es::Matrix{T}, + kernel::RK + ) where {T <: AbstractFloat, RK <: ReproducingKernel_1} + cond = T(0.0) + gram = _get_gram(nodes, kernel) + try + evs = svdvals!(gram) + maxevs = maximum(evs) + minevs = minimum(evs) + if minevs > T(0.0) + cond = maxevs / minevs + cond = 10.0^floor(log10(cond)) + end + catch + end + return cond +end + +# ``` +# Get estimation of the Gram matrix condition number +# Brás, C.P., Hager, W.W. & Júdice, J.J. An investigation of feasible descent algorithms for estimating the condition number of a matrix. TOP 20, 791–809 (2012). +# https://link.springer.com/article/10.1007/s11750-010-0161-9 +# ``` +function _estimate_cond(gram::Matrix{T}, + chol::LinearAlgebra.Cholesky{T,Array{T,2}}, + nit = 3 + ) where T <: AbstractFloat + if isnothing(gram) + throw(DomainError(gram, "Parameter `gram` is `nothing'.")) + end + if isnothing(chol) + throw(DomainError(chol, "Parameter `chol` is `nothing'.")) + end + mat_norm = norm(gram, 1) + n = size(gram, 1) + x = Vector{T}(undef, n) + @. x = T(1.0) / T(n) + z = Vector{T}(undef, n) + gamma = T(0.0) + for it = 1:nit + z = ldiv!(z, chol, x) + gamma = T(0.0); + for i = 1:n + gamma += abs(z[i]) + z[i] = sign(z[i]) + end + z = ldiv!(z, chol, copy(z)) + zx = z ⋅ x + idx = 1 + for i = 1:n + z[i] = abs(z[i]) + if z[i] > z[idx] + idx = i + end + end + if z[idx] <= zx + break + end + @. x = T(0.0) + x[idx] = T(1.0) + end + cond = T(10.0)^floor(log10(mat_norm * gamma)) + return cond +end