Revise logarithmic_norm code and tests

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 
 [compat]
 Documenter = "1"
+DocumenterCitations = "1.3"

docs/make.jl

@@ -1,4 +1,6 @@
-using Documenter, ReachabilityBase
+using Documenter, ReachabilityBase, DocumenterCitations
+
+bib = CitationBibliography(joinpath(@__DIR__, "src", "refs.bib"); style=:alpha)
 
 DocMeta.setdocmeta!(ReachabilityBase, :DocTestSetup,
                     :(using ReachabilityBase); recursive=true)
@@ -12,6 +14,7 @@ makedocs(; sitename="ReachabilityBase.jl",
                   ReachabilityBase.Distribution, ReachabilityBase.Subtypes,
                   ReachabilityBase.Arrays, ReachabilityBase.Timing],
          pagesonly=true,
+         plugins=[bib],
          pages=["Home" => "index.md",
                 "Library" => Any[
                                  #
@@ -28,6 +31,7 @@ makedocs(; sitename="ReachabilityBase.jl",
                                  "Basetype" => "lib/Basetype.md"
                                  #
                                  ],
+                "Bibliography" => "bibliography.md",
                 "About" => "about.md"])
 
 deploydocs(; repo="github.com/JuliaReach/ReachabilityBase.jl.git",

docs/src/bibliography.md

@@ -0,0 +1,4 @@
+# Bibliography
+
+```@bibliography
+```

docs/src/lib/Arrays.md

@@ -29,6 +29,7 @@ ismultiple
 ispermutation
 issquare
 iswellconditioned
+logarithmic_norm
 matrix_type
 nonzero_columns
 nonzero_indices

docs/src/refs.bib

@@ -0,0 +1,6 @@
+@book{DesoerV09,
+  title={Feedback Systems: Input-Output Properties},
+  author={Desoer, Charles A. and Vidyasagar, Mathukumalli},
+  year={2009},
+  publisher={SIAM}
+}

src/Arrays/Arrays.jl

@@ -6,7 +6,7 @@ This module provides machinery for vectors and matrices.
 module Arrays
 
 using Base: \, /, transpose
-using LinearAlgebra: Diagonal, cond, det, diag, eigvals, norm, qr
+using LinearAlgebra: Diagonal, Hermitian, cond, det, diag, eigmax, norm, qr
 using Random: AbstractRNG, GLOBAL_RNG
 using Requires: @require
 using SparseArrays: AbstractSparseArray, AbstractSparseMatrix,

src/Arrays/logarithmic_norm.jl

@@ -1,27 +1,38 @@
-# =======================================================================
-# Logarithmic norms of matrices (a.k.a. matrix measures)
-#
-# Reference:
-#
-# C. Desoer and M. Vidyasagar, Feedback Systems: Input-Output Properties.
-# Society for Industrial and Applied Mathematics, 2009.
-# =======================================================================
-
-# logarithmic norm (also known as matrix measure) for commonly used p-norms
+"""
+    logarithmic_norm(A::AbstractMatrix, p::Real=Inf)
+
+Compute the logarithmic norm (also known as matrix measure) for commonly used p-norms.
+
+### Input
+
+- `A` -- matrix
+- `p` -- (optional, default: `Inf`) p-norm
+
+### Output
+
+A number representing the logarithmic norm.
+
+### Notes
+
+The implementation currently supports the following values for `p`: `1`, `2`, `Inf`.
+
+See [DesoerV09](@citet) for a reference.
+"""
 function logarithmic_norm(A::AbstractMatrix, p::Real=Inf)
     if p == Inf
-        return logarithmic_norm_inf(A)
+        return _logarithmic_norm_inf(A)
     elseif p == 1
-        return logarithmic_norm_1(A)
+        return _logarithmic_norm_1(A)
     elseif p == 2
-        return logarithmic_norm_2(A)
+        return _logarithmic_norm_2(A)
     else
-        throw(ArgumentError("logarithmic norm only implemented for p = 1, 2, or Inf, got p = $p"))
+        throw(ArgumentError("`logarithmic_norm` is only implemented for p = 1, 2, or Inf, " *
+                            "but got p = $p"))
     end
 end
 
 # max_j a_jj + ∑_{i ≠ j} |a_ij|
-function logarithmic_norm_1(A::AbstractMatrix{N}) where {N}
+function _logarithmic_norm_1(A::AbstractMatrix)
     out = -Inf
     @inbounds for j in axes(A, 2)
         α = A[j, j]
@@ -38,7 +49,7 @@ function logarithmic_norm_1(A::AbstractMatrix{N}) where {N}
 end
 
 # max_i a_ii + ∑_{j ≠ i} |a_ij|
-function logarithmic_norm_inf(A::AbstractMatrix{N}) where {N}
+function _logarithmic_norm_inf(A::AbstractMatrix)
     out = -Inf
     @inbounds for i in axes(A, 1)
         α = A[i, i]
@@ -54,9 +65,9 @@ function logarithmic_norm_inf(A::AbstractMatrix{N}) where {N}
     return out
 end
 
-# max_j  1/2 * λⱼ(A + A^T)
-function logarithmic_norm_2(A::AbstractMatrix)
-    B = A + A'
-    λ = eigvals(B)
-    return maximum(λ) / 2
+# max_j  1/2 * λⱼ(A + Aᵀ)
+function _logarithmic_norm_2(A::AbstractMatrix)
+    B = Hermitian(A + A')  # type information that all eigenvalues are real
+    λ = eigmax(B)
+    return λ / 2
 end

test/Arrays/logarithmic_norm.jl

@@ -1,23 +1,26 @@
-using ReachabilityBase.Arrays
+using ReachabilityBase.Arrays, Test
+using ReachabilityBase.Comparison: isidentical
 
-@testset "Logarithmic norm" begin
-    # Example 1
-    A = [1 2; 3 -4]
-
-    @test logarithmic_norm(A, 1) == 4
-    @test logarithmic_norm(A, Inf) == 3
-
-    # Example 2 (see Example 2 in "Explicit Error Bounds for Carleman Linearization.
-    # Marcelo Forets and Amaury Pouly (2018)").
-    A = [0 1; -1 -2]
+for N in [Float64, Float32, Rational{Int}]
+    A = N[1 2; 3 -4]
+    @test_throws ArgumentError logarithmic_norm(A, 3)
+    @test isidentical(logarithmic_norm(A), N(3))  # default norm is `Inf`
+    @test isidentical(logarithmic_norm(A, Inf), N(3))
+    @test isidentical(logarithmic_norm(A, 1), N(4))
+    # see Example 2 in:
+    # "Marcelo Forets and Amaury Pouly: Explicit Error Bounds for Carleman Linearization (2018)"
+    A = N[0 1; -1 -2]
     Aop2 = kron(A, A)
+    @test isidentical(logarithmic_norm(A, Inf), N(1))
+    @test isidentical(logarithmic_norm(Aop2, Inf), N(9))
+    @test isidentical(logarithmic_norm(A, 1), N(1))
+    @test isidentical(logarithmic_norm(Aop2, 1), N(9))
+end
 
-    @test logarithmic_norm(A, Inf) == 1
-    @test logarithmic_norm(Aop2, Inf) == 9
-
-    @test logarithmic_norm(A, 1) == 1
-    @test logarithmic_norm(Aop2, 1) == 9
-
-    @test logarithmic_norm(A, 2) == 0
-    @test logarithmic_norm(Aop2, 2) ≈ 4.23606797749979
+for N in [Float64, Float32]
+    A = N[0 1; -1 -2]  # see above
+    Aop2 = kron(A, A)
+    @test isidentical(logarithmic_norm(A, 2), N(0))
+    res = logarithmic_norm(Aop2, 2)
+    @test res isa N && res ≈ 4.23606797749979
 end

test/runtests.jl

@@ -11,9 +11,6 @@ seed!(1234)  # fix RNG seed for reproducibility
     @testset "arrays" begin
         include("Arrays/arrays.jl")
     end
-    @testset "logarithmic_norm" begin
-        include("Arrays/logarithmic_norm.jl")
-    end
     @testset "matrix_operations" begin
         include("Arrays/matrix_operations.jl")
     end
@@ -26,6 +23,9 @@ seed!(1234)  # fix RNG seed for reproducibility
     @testset "SingleEntryVector" begin
         include("Arrays/SingleEntryVector.jl")
     end
+    @testset "logarithmic_norm" begin
+        include("Arrays/logarithmic_norm.jl")
+    end
 end
 @testset "Comparison" begin
     include("Comparison/comparison.jl")