export, reformat, add docstrings

Author: dpo

src/DiagonalHessianApproximation.jl

@@ -1,3 +1,5 @@
+export DiagonalQN, SpectralGradient
+
 function mulSquareOpDiagonal!(res, d, v, α, β::T) where {T <: Real}
   if β == zero(T)
     res .= α .* d .* v
@@ -13,7 +15,8 @@ Andrei, N.
 A diagonal quasi-Newton updating method for unconstrained optimization. 
 https://doi.org/10.1007/s11075-018-0562-7
 """
-mutable struct DiagonalQN{T <: Real, I <: Integer, V <: AbstractVector{T}} <: AbstractDiagonalQuasiNewtonOperator{T} 
+mutable struct DiagonalQN{T <: Real, I <: Integer, V <: AbstractVector{T}} <:
+               AbstractDiagonalQuasiNewtonOperator{T}
   d::V # Diagonal of the operator
   nrow::I
   ncol::I
@@ -30,43 +33,54 @@ mutable struct DiagonalQN{T <: Real, I <: Integer, V <: AbstractVector{T}} <: Ab
   allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
-DiagonalQN(d::AbstractVector{T}) where {T <: Real} = 
-  DiagonalQN(
-    d,
-    length(d),
-    length(d),
-    true, 
-    true,  
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
-    0,
-    0,
-    0,
-    true,
-    true,
-    true)
+"""
+    DiagonalQN(d)
+
+Construct a linear operator that represents a diagonal quasi-Newton approximation.
+The approximation satisfies the weak secant equation and is not guaranteed to be
+positive definite.
+
+# Arguments
+
+- `d::AbstractVector`: initial diagonal approximation.
+"""
+DiagonalQN(d::AbstractVector{T}) where {T <: Real} = DiagonalQN(
+  d,
+  length(d),
+  length(d),
+  true,
+  true,
+  (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
+  (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
+  (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
+  0,
+  0,
+  0,
+  true,
+  true,
+  true,
+)
 
 # update function
 # s = x_{k+1} - x_k
 # y = ∇f(x_{k+1}) - ∇f(x_k)
 function push!(
-  B::DiagonalQN{T,I,V},
+  B::DiagonalQN{T, I, V},
   s::V,
-  y::V
-  ) where {T <: Real, I <: Integer, V <: AbstractVector{T}}
+  y::V,
+) where {T <: Real, I <: Integer, V <: AbstractVector{T}}
   trA2 = zero(T)
   for i in eachindex(s)
     trA2 += s[i]^4
   end
-  sT_s = dot(s,s)
-  sT_y = dot(s,y)
+  sT_s = dot(s, s)
+  sT_y = dot(s, y)
   sT_B_s = sum(s[i]^2 * B.d[i] for i ∈ eachindex(s))
   if trA2 == 0
     error("Cannot divide by zero and trA2 = 0")
   end
-  q = (sT_y + sT_s - sT_B_s)/trA2
-  B.d .+= q .* s.^2 .- 1
+  q = (sT_y + sT_s - sT_B_s) / trA2
+  B.d .+= q .* s .^ 2 .- 1
   return B
 end
 
@@ -77,7 +91,7 @@ Birgin, E. G., Martínez, J. M., & Raydan, M.
 Spectral Projected Gradient Methods: Review and Perspectives. 
 https://doi.org/10.18637/jss.v060.i03
 """
-mutable struct SpectralGradient{T <: Real, I <: Integer} <: AbstractDiagonalQuasiNewtonOperator{T} 
+mutable struct SpectralGradient{T <: Real, I <: Integer} <: AbstractDiagonalQuasiNewtonOperator{T}
   d::T # Diagonal coefficient of the operator (multiple of the identity)
   nrow::I
   ncol::I
@@ -94,34 +108,48 @@ mutable struct SpectralGradient{T <: Real, I <: Integer} <: AbstractDiagonalQuas
   allocated5::Bool # true for 5-args mul!, false for 3-args mul! until the vectors are allocated
 end
 
-SpectralGradient(d::T, n::I) where {T <: Real, I <: Integer} = 
+"""
+        SpectralGradient(σ, n)
+
+Construct a spectral gradient Hessian approximation.
+The approximation is defined as σI.
+
+# Arguments
+
+- `σ::Real`: initial positive multiple of the identity;
+- `n::Int`: operator size.
+"""
+function SpectralGradient(d::T, n::I) where {T <: Real, I <: Integer}
+  @assert d > 0
   SpectralGradient(
     d,
     n,
     n,
     true,
-    true, 
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
-    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β), 
+    true,
+    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
+    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
+    (res, v, α, β) -> mulSquareOpDiagonal!(res, d, v, α, β),
     0,
     0,
     0,
     true,
     true,
-    true)
+    true,
+  )
+end
 
 # update function
 # s = x_{k+1} - x_k
 # y = ∇f(x_{k+1}) - ∇f(x_k)
 function push!(
-  B::SpectralGradient{T,I},
+  B::SpectralGradient{T, I},
   s::V,
-  y::V
-  ) where {T <: Real, I <: Integer, V <: AbstractVector{T}}
+  y::V,
+) where {T <: Real, I <: Integer, V <: AbstractVector{T}}
   if all(s .== 0)
     error("Cannot divide by zero and s .= 0")
   end
-  B.d = dot(s,y)/dot(s,s)
+  B.d = dot(s, y) / dot(s, s)
   return B
-end
+end