Keeping well-defined mathematical functions orthogonal to a project structs design

General Usage
,
1

I have a doubt regarding a design pattern. To expose it I’m going to refer to IterativeSolvers.jl’s Conjugate Gradient algorithm.

For the end user, they supply the function cg(A, b; kwargs...), which solves the equation Ax = b and returns the solution x. To solve it, the package internally defines a mutable struct containing all vectors and matrices that will be needed in the iterative algorithm,

mutable struct CGIterable{matT, solT, vecT, numT <: Real}
    A::matT
    x::solT
    r::vecT
    c::vecT
    u::vecT
    tol::numT
    residual::numT
    prev_residual::numT
    maxiter::Int
    mv_products::Int
end

and essentially all the work is done in the following function, which performs a step of the conjugate gradient algorithm, updating the struct of type CGIterable:

###############
# Ordinary CG #
###############

function iterate(it::CGIterable, iteration::Int=start(it))
    # Check for termination first
    if done(it, iteration)
        return nothing
    end

    # u := r + βu (almost an axpy)
    β = it.residual^2 / it.prev_residual^2
    it.u .= it.r .+ β .* it.u

    # c = A * u
    mul!(it.c, it.A, it.u)
    α = it.residual^2 / dot(it.u, it.c)

    # Improve solution and residual
    it.x .+= α .* it.u
    it.r .-= α .* it.c

    it.prev_residual = it.residual
    it.residual = norm(it.r)

    # Return the residual at item and iteration number as state
    it.residual, iteration + 1
end

The problem I see with this is that CG is a mathematically well-defined algorithm which needs a specific set of matrices and vectors as input to produce an output. However, the implementation above depends on the particular design of the struct of type CGIterable, and can only be used if this object is built.

For me it seems more natural to have instead a function CG_step that is independent of the particular design of the structs in my project, such as:

function CG_step(x, A, r, c, u, residual, prev_residual)
    # Check for termination first
    # if done(it, iteration)
    #     return nothing
    # end

    # u := r + βu (almost an axpy)
    β = residual^2 / prev_residual^2
    u .= r .+ β .* u

    # c = A * u
    mul!(c, A, u)
    α = residual^2 / dot(u, c)

    # Improve solution and residual
    x .+= α .* u
    r .-= α .* c

    prev_residual = residual
    residual = norm(r)

    # Return the residual at item and iteration number as state
    residual, iteration + 1
end

and then have an additional function so that I can use my structs throughout my project:

function iterate(it::CGIterable, iteration::Int=start(it))
    return CG_step(it.x, it.A, it.r, it.c, it.u, it.residual, it.prev_residual)
end

That is, what I am trying to convey is that I would think it is desirable to make well-defined mathematical functions that do all the job orthogonal to the particular design of project’s structs. This way, if I ever want to change my structs, I don’t need to touch the actually important function (which could potentially introduce bugs), but just the interface function defined just above. The actual mathematical function would be coded “once and for all” in Julia.

Of course, I guess that I am missing some important points that make this approach undesirable, but I wanted to ask about your opinions on this. What I said makes any sense? What are the advantages/disadvantages of doing what I said? Why didn’t they code the algorithm this way?

4 Likes
2

Your version can work well in some cases, but look at the downsides:

  • You need to pass a large number of parameters to your CG_step. This could get unwieldy in other cases.
  • How can the function store the changes it makes to the parameters? In-place assignments of arrays will work, but updating scalars will not work. You can play with multiple return values as you do with residual (but you don’t use it actually?). But this makes things more brittle and complicated.

What’s best depends on the particular problem, but in general I think wrapping many related values in a struct is good software engineering practice.

And tangentially related is this quote from Linus Torvalds (see discussion here):

Bad programmers worry about the code. Good programmers worry about data structures and their relationships.

It’s often a good idea to prioritize the choice of data structure even if doing so you lose some mathematical purity.

2 Likes
3

I guess it’s a good point. My biggest concern was that, if I wanted to check my own implementation of the algorithm with IterativeSolvers.jl’s one, I would need to first create their data structures to be able to call the function, so that what I gained in passing fewer inputs is lost because I still need to create a struct beforehand.

However, I agree that the code is much cleaner from the developer point of view.

1 Like
4

Multiple dispatch can help with this. While you could write different versions of the function for different structs, you could also write adapters.

function CG_step(x, A, r, c, u, residual, prev_residual; maxiter, mv_products)
   cgiter = CGIterable(A, x, r, c, u, tol, residual, prev_residual, maxiter, mv_products)
   return IterativeSolvers.iterate(cgiter)
end

function IterativeSolvers.iterate(it::MyCGIterable)
    cgiter = convert(IterativeSolvers.CGIterable, it)
    return IterativeSolvers.iterate(cgiter)
end

function Base.convert(::Type{IterativeSolvers.CGIterable}, it::MyCGIterable)
    return IterativeSolvers.CGIterable( ... )
end
1 Like