rank estimation clean up, add circles method

Author: njericha

src/BlockTensorFactorization.jl

@@ -83,14 +83,14 @@ export getconstraint, smart_insert!, smart_interlace!, group_by_factor
 
 # High level / user-interface
 #include("./factorize.jl")
-export factorize
+export factorize, default_kwargs
 
 #include("./multiscale.jl")
 export multiscale_factorize
 export coarsen, interpolate, scale_constraint
 
 #include("./rankdetection.jl")
 export rank_detect_factorize
-export possible_ranks
+export max_possible_rank
 
 end # module BlockTensorFactorization

src/Core/Core.jl

@@ -93,14 +93,14 @@ export getconstraint, smart_insert!, smart_interlace!, group_by_factor
 
 # High level / user-interface
 include("./factorize.jl")
-export factorize
+export factorize, default_kwargs
 
 include("./multiscale.jl")
 export multiscale_factorize
 export coarsen, interpolate, scale_constraint
 
 include("./rankdetection.jl")
 export rank_detect_factorize
-export possible_ranks
+export max_possible_rank
 
 end # module Core

src/Core/curvaturetools.jl

@@ -296,35 +296,37 @@ end
 
 """Extracts the first and second derivatives of the splines at the knots"""
 function d_dx_and_d2_dx2_spline(y::AbstractVector{<:Real}; h=1)
-    _, b, c, _ = cubic_spline_coefficients(y::AbstractVector{<:Real}; h=1)
+    _, b, c, _ = cubic_spline_coefficients(y::AbstractVector{<:Real}; h)
     dy_dx = c
     dy2_dx2 = 2b
     return dy_dx, dy2_dx2
 end
 
 
 """
-    curvature(y::AbstractVector{<:Real})
+    curvature(y::AbstractVector{<:Real}; method=:finite_differences)
 
 Approximates the signed curvature of a function given evenly spaced samples.
 
 Uses [`d_dx`](@ref) and [`d2_dx2`](@ref) to approximate the first two derivatives.
 """
 function curvature(y::AbstractVector{<:Real}; method=:finite_differences, kwargs...)
-    if method == finite_differences
+    if method == :finite_differences
         dy_dx = d_dx(y; kwargs...)
         dy2_dx2 = d2_dx2(y; kwargs...)
         return @. dy2_dx2 / (1 + dy_dx^2)^1.5
     elseif method == :splines
         dy_dx, dy2_dx2 = d_dx_and_d2_dx2_spline(y; h=1)
         return @. dy2_dx2 / (1 + dy_dx^2)^1.5
+    elseif method == :circles
+        return circumscribed_standard_curvature(y)
     else
         throw(ArgumentError("method $method not implemented"))
     end
 end
 
 """
-    standard_curvature(y::AbstractVector{<:Real})
+    standard_curvature(y::AbstractVector{<:Real}; method=:finite_differences)
 
 Approximates the signed curvature of a function, scaled to the unit box ``[0,1]^2``.
 
@@ -341,32 +343,62 @@ function standard_curvature(y::AbstractVector{<:Real}; method=:finite_difference
         # y_max = 1
         dy_dx, dy2_dx2 = d_dx_and_d2_dx2_spline(y; h=Δx)
         return @. dy2_dx2 / (1 + dy_dx^2)^1.5
+    elseif method == :circles
+        return circumscribed_standard_curvature(y)
     else
         throw(ArgumentError("method $method not implemented"))
     end
 end
 
-"""
-Finds the radius of the circumscribed circle between points (a,f), (b,g), (c,h)
-"""
-function circumscribed_radius((a,f),(b,g),(c,h))
-    d = 2*(a*(g-h)+b*(h-f)+c*(f-g))
-    p = ((a^2+f^2)*(g-h)+(b^2+g^2)*(h-f)+(c^2+h^2)*(f-g)) / d
-    q = ((a^2+f^2)*(b-c)+(b^2+g^2)*(c-a)+(c^2+h^2)*(a-b)) / d
-    r = sqrt((a-p)^2+(f-q)^2)
-    return r
-end
+# """
+# Finds the radius of the circumscribed circle between points (a,f), (b,g), (c,h)
+# """
+# function circumscribed_radius((a,f),(b,g),(c,h))
+#     d = 2*(a*(g-h)+b*(h-f)+c*(f-g))
+#     p = ((a^2+f^2)*(g-h)+(b^2+g^2)*(h-f)+(c^2+h^2)*(f-g)) / d
+#     q = ((a^2+f^2)*(b-c)+(b^2+g^2)*(c-a)+(c^2+h^2)*(a-b)) / d
+#     r = sqrt((a-p)^2+(f-q)^2)
+#     return r
+# end
 
 function circumscribed_standard_curvature(y)
-    n = length(v)
+    n = length(y)
     ymax = maximum(y)
     y = y / ymax
-    k = zero(ymax)
+    k = zero(y)
     a, b, c = 0, 1/n, 2/n
     for i in eachindex(k)[2:end-1]
-        k[i] = 1 / circumscribed_radius((a,y[i-1]),(b,y[i]),(c,y[i+1]))
+        k[i] = signed_circle_curvature((a,y[i-1]),(b,y[i]),(c,y[i+1]))
+        #k[i] = 1 / circumscribed_radius((a,y[i-1]),(b,y[i]),(c,y[i+1]))
     end
     k[1] = k[2]
     k[end] = k[end-1]
     return k
 end
+
+"""radius r and center point (p,q) of the circle passing through the three points"""
+function three_point_circle((a,f),(b,g),(c,h))
+    fg = f-g
+    gh = g-h
+    hf = h-f
+    ab = a-b
+    bc = b-c
+    ca = c-a
+    a2 = a^2
+    b2 = b^2
+    c2 = c^2
+    f2 = f^2
+    g2 = g^2
+    h2 = h^2
+    p = (a2*gh + b2*hf + c2*fg - gh*hf*fg) / (a*gh + b*hf + c*fg) / 2
+    q = (f2*bc + g2*ca + h2*ab - bc*ca*ab) / (f*bc + g*ca + h*ab) / 2
+    r = sqrt((a-p)^2 + (f-q)^2)
+    return r, (p, q)
+end
+
+function signed_circle_curvature((a,f),(b,g),(c,h))
+    @assert a < b < c
+    r, _ = three_point_circle((a,f),(b,g),(c,h))
+    sign = g > (f+h)/2 ? -1 : 1
+    return sign / r
+end

src/Core/rankdetection.jl

@@ -1,11 +1,26 @@
 """
-    rank_detect_factorize(Y; online_rank_estimation=false, rank=nothing, model=Tucker1, kwargs...)
+    rank_detect_factorize(Y; kwargs...)
 
 Wraps `factorize()` with rank detection.
 
 Selects the rank that maximizes the standard curvature of the Relative Error (as a function of rank).
+
+# Keywords
+- `online_rank_estimation`: `false`. Set to `true` to stop testing larger ranks after the first peak in curvature
+- `curvature_method`: `:splines`. Can also pick `:finite_differences` (faster but less accurate) or `circles` (fastest and smallest memory but more sensitive to results from `factorize`)
+- `model`: `Tucker1`. Only rank detection with `Tucker1` and `CPDecomposition` is currently implemented
+- `max_rank`: `max_possible_rank(Y, model)`. Test ranks from `1` up to `max_rank`. Defaults to largest possible rank under the model
+- `rank`: `nothing`. If a rank is passed, rank detection is ignored and `factorize(Y; kwargs...)` is called
+
+Any other keywords from [`factorize`](@ref), full list given by [`default_kwargs`](@ref).
 """
-function rank_detect_factorize(Y; online_rank_estimation=false, rank=nothing, model=Tucker1, kwargs...)
+function rank_detect_factorize(Y;
+    online_rank_estimation=false,
+    rank=nothing,
+    model=Tucker1,
+    max_rank=max_possible_rank(Y, model),
+    curvature_method=:splines,
+    kwargs...)
     if isnothing(rank)
         # Initialize output and final error lists
         all_outputs = []
@@ -20,8 +35,10 @@ function rank_detect_factorize(Y; online_rank_estimation=false, rank=nothing, mo
             kwargs[:stats] = [RelativeError, kwargs[:stats]...] # not using pushfirst! since kwargs[:stats] could be a Tuple
         end
         kwargs[:model] = model # add the model back into kwargs
+        kwargs[:curvature_method] = curvature_method# :splines
+        kwargs[:max_rank] = max_rank
 
-        for rank in possible_ranks(Y, model)
+        for rank in 1:max_rank
             @info "Trying rank=$rank..."
 
             kwargs[:rank] = rank # add the rank into kwargs
@@ -35,7 +52,7 @@ function rank_detect_factorize(Y; online_rank_estimation=false, rank=nothing, mo
             @info "Final relative error = $final_rel_error"
 
             if (online_rank_estimation == true) && length(final_rel_errors) >= 3 # Need at least 3 points to evaluate curvature
-                curvatures = standard_curvature(final_rel_errors; method=:finite_differences) # method=:splines
+                curvatures = standard_curvature(final_rel_errors; method=curvature_method) # method=:splines
                 if curvatures[end] ≈ maximum(curvatures) # want the last curvature to be significantly smaller than the max
                     continue
                 else
@@ -57,9 +74,9 @@ function rank_detect_factorize(Y; online_rank_estimation=false, rank=nothing, mo
 end
 
 """
-    possible_ranks(Y, model)
+    max_possible_rank(Y, model)
 
-Returns the rank of possible ranks `Y` could have under the `model`.
+Returns the maximum rank possible `Y` could have under the `model`.
 
 For matrices `I × J` this is `1:min(I, J)`. This is can be extended to tensors for different type
 of decompositions.
@@ -70,16 +87,16 @@ The CP-rank is `≤ minimum_{n} (prod(I1,...,IN) / In)` for tensors `I1 × … 
 some shapes have have tighter upper bounds. For example, `2 × I × I` tensors over ℝ have a maximum
 rank of `floor(3I/2)`.
 """
-function possible_ranks(Y, model)
+function max_possible_rank(Y, model) # TODO store this info with the Corresponding AbstractDecomposition
     if model <: Tucker1
         I, Js... = size(Y)
         max_rank = min(I, prod(Js))
-        return 1:max_rank
+        return max_rank
     elseif model <: CPDecomposition
         Is = size(Y)
         # There exist tighter upper bounds for particular shapes like I×I×K, but this a simple upper bound that works for all shapes
         max_rank = minimum(prod(Is) .÷ Is) # ÷ is Integer division
-        return 1:max_rank
+        return max_rank
     else
         error("Possible ranks for models of type $model are not implemented")
     end

test/Project.toml

@@ -2,4 +2,5 @@
 BlockTensorFactorization = "07b766a1-0096-4e53-b687-1cd07becaccb"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/runtests.jl

@@ -7,6 +7,7 @@ using Test
 using Random
 Random.seed!(3141592653589) # Reproducibility of random initializations
 using LinearAlgebra
+using Statistics
 
 using BlockTensorFactorization
 
@@ -803,10 +804,39 @@ end
 end
 
 @testset "RankDetection" begin
+    @test BlockTensorFactorization.Core.three_point_circle((1,2), (2,1), (5,2)) == (√5, (3, 3))
+
+    # any smooth function s.t.
+    # f(0) = 1, f(1) = 0, f''(1) = 0
+    f(x) = 3x^3 - 5x^2 + x + 1
+    # k(x) = f''(x) / (1 + (f'(x))^2)^(3/2)
+    k(x) = (18x-10)/(1 + (1 - 10x + 9x^2)^2)^(1.5) # closed form true curvature
+    x = range(0, 1, length=20)[2:end] # exclude x=0 point (necessary for :splines)
+    y = f.(x)
+    k_true = k.(x)
+    methods = (:finite_differences, :splines, :circles)
+    k_finite_differences, k_splines, k_circles = (standard_curvature(y; method=method) for method in methods)
+    # Mean Absolute Percentage Error
+    MAPE(test_vals, true_vals) = mean(@. abs((test_vals - true_vals)/true_vals))
+    @test MAPE(k_splines, k_true) < 0.07 # 7%
+    @test MAPE(k_circles, k_true) < 0.08 # 8%
+    @test MAPE(k_finite_differences, k_true) < 0.09 # 9%
+
     T = Tucker1((10, 10, 10), 3)
     Y = array(T)
-    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1)
+    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:splines)
+    @test kwargs[:rank] == 3
+
+    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:circles)
     @test kwargs[:rank] == 3
+
+    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=Tucker1, curvature_method=:finite_differences)
+    @test kwargs[:rank] == 3
+
+    T = CPDecomposition((10, 11, 12), 4)
+    Y = array(T)
+    decomposition, stats, kwargs, final_rel_errors = rank_detect_factorize(Y; model=CPDecomposition, curvature_method=:splines, online_rank_estimation=true)
+    @test kwargs[:rank] == 4
 end
 
 end