PartitionedStructures.jl: Tutorial

This tutorial shows how PartitionedStructures.jl can define partitioned quasi-Newton approximations, by using partitioned vectors and partitioned matrices. Those partitioned structures are strongly related to partially-separable functions.

What are the partially-separable structure and the partitioned quasi-Newton approximations

The partitioned quasi-Newton methods exploit the partially-separable structure of $f:\R^n \to \R$

\[ f(x) = \sum_{i=1}^N f_i (U_i x) , \; f_i : \R^{n_i} \to \R, \; U_i \in \R^{n_i \times n},\; n_i \ll n,\]

as the sum of the element functions $f_i$. The gradient $\nabla f$ and the Hessian $\nabla^2 f$

\[\nabla f(x) = \sum_{i=1}^N U_i^\top f_i (U_i x), \quad \nabla^2 f(x) = \sum_{i=1}^N U_i^\top f_i (U_i x) U_i,\]

accumulate the element derivatives $\nabla f_i$ and $\nabla^2 f_i$.

The partitioned structure of the Hessian allows the definition of partitioned quasi-Newton approximations $B \approx \nabla^2 f$, such that $B$ accumulates every element Hessian approximation $B_i \approx \nabla^2 f_i$

\[B = \sum_{i=1}^N U_i^\top B_i U_i\]

The partitioned quasi-Newton approximations respect the sparsity structure of $\nabla^2 f$, which is not the case of classical quasi-Newton approximation (e.g. BFGS, SR1). Moreover, the rank of the partitioned updates may be proportional to the number of elements $N$, whereas classical quasi-Newton approximations are low-rank updates. A partitioned quasi-Newton update may update every element-Hessian $B_i$ at each step $s$. It requires $B_i$, $U_i s$ and $\nabla f_i (U_i (x+s)) - \nabla f_i (U_i x)$, and therefore we have to store such an approximation and vectors for every element.

Reference

• A. Griewank and Ph. L. Toint. Partitioned variable metric updates for large structured optimization problems. Numer. Math., 39:119–137, 1982, 10.1007/BF01399316

Example: the partitioned structure of a quadratic

Let's take the quadratic f as an example

f(x) = x[1]^2 + x[1]*x[2] + x[2]^2 + x[3]^2 + 3x[2]*x[3]

f (generic function with 1 method)

f can be considered as the sum of two element functions

f1(y) = y[1]^2 + y[1]*y[2]
f2(y) = y[1]^2 + y[2]^2 + 3y[1]*y[2]

f2 (generic function with 1 method)

Define $U_1$ and $U_2$ indicating which variables are required by each element function:

U1 = [1 0 0; 0 1 0]
U2 = [0 1 0; 0 0 1]

2×3 Matrix{Int64}:
 0  1  0
 0  0  1

However, dense matrices produce a memory issue for large problems. Instead, we use a linear operator only informing the indices of the variables

U1 = [1, 2]
U2 = [2, 3]

2-element Vector{Int64}:
 2
 3

By gathering the different $U_i$ together

U = [U1, U2]

2-element Vector{Vector{Int64}}:
 [1, 2]
 [2, 3]

we define the function exploiting the partially-separable structure f_pss as

f_pss(x, U) = f1(x[U[1]]) + f2(x[U[2]])
x0 = [2., 3., 4.]

f(x0) == f_pss(x0, U)

true

Similarly, we can compute: the gradient, the element gradients and explicit how the gradient is partitioned

∇f(x) = [2x[1] + x[2], x[1] + 2x[2] + 3x[3], 2x[3] + 3x[2]]
∇f1(y) = [2y[1] + y[2], y[1]]
∇f2(y) = [2y[1] + 3y[2], 2y[2] + 3y[1]]
function ∇f_pss(x, U)
  gradient = zeros(length(x))
  gradient[U1] .+= ∇f1(x[U[1]])
  gradient[U2] .+= ∇f2(x[U[2]])
  return gradient
end
∇f(x0) == ∇f_pss(x0, U)

true

However, ∇f_pss accumulates directly the element gradients and does not store the value of each element gradients ∇f1, ∇f2. We would like to store every element gradient, so we can later update partitioned quasi-Newton approximations. Thus, we define a partitioned vector from U to store each element gradient and form $\nabla f$ when required:

using PartitionedStructures
U = [U1, U2]
n = length(x0)
# create the partitioned vector
partitioned_gradient_x0 = create_epv(U)

Elemental_pv{Float64}(2, 3, Elemental_elt_vec{Float64}[Elemental_elt_vec{Float64}([0.26530257356243503, 0.17482000587233903], [1, 2], 2), Elemental_elt_vec{Float64}([0.9611384261051259, 0.8918990543826427], [2, 3], 2)], [0.0, 0.0, 0.0], [[1], [1, 2], [2]], [1, 2, 3])

We set the value of each element vector to the corresponding element gradient

# return every element gradient
vector_gradient_element(x, U) = [∇f1(x[U[1]]), ∇f2(x[U[2]])] :: Vector{Vector{Float64}}
# set each element vector to its corresponding element gradient
set_epv!(partitioned_gradient_x0, vector_gradient_element(x0, U))

# Build the gradient vector
build_v!(partitioned_gradient_x0)

# with the same value as the gradient
get_v(partitioned_gradient_x0) == ∇f(x0)

true

Approximate the Hessian $\nabla^2 f$

There are at least two ways to approximate $\nabla^2 f$: a classical quasi-Newton approximation (ex: BFGS) and a partitioned quasi-Newton approximation (ex: PBFGS). This section presents both methods and exposes the sparse structure of the partitioned variant.

Quasi-Newton approximation of the quadratic (BFGS)

In the case of the BFGS method, we want to approximate the Hessian from s = x1 - x0

x1 = [1., 2., 2.]
s = x1 .- x0

3-element Vector{Float64}:
 -1.0
 -1.0
 -2.0

the gradient difference y

y = ∇f(x1) .- ∇f(x0)

3-element Vector{Float64}:
 -3.0
 -9.0
 -7.0

and the approximation B, initially set to the identity

B = [ i==j ? 1. : 0. for i in 1:n, j in 1:n]

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

By applying the BFGS update, you satisfy the secant equation Bs = y

# PartitionedStructures.jl exports a BFGS update
B_BFGS = BFGS(s,y,B)

using LinearAlgebra
# numerical verification of the secant equation
atol = sqrt(eps(eltype(s)))
norm(B_BFGS * s - y) < atol

true

but the approximation B_BFGS is dense:

B_BFGS

3×3 Matrix{Float64}:
 1.17949   0.871795  0.474359
 0.871795  3.94872   2.08974
 0.474359  2.08974   2.21795

which makes it impractical for large-scale problems.

Partitioned quasi-Newton approximation of the quadratic (PBFGS)

In order to make a sparse quasi-Newton approximation of $\nabla^2 f$, you may define a partitioned matrix with the same partially-separable structure as partitioned_gradient_x0 where each element matrix is set to the identity

partitioned_matrix = epm_from_epv(partitioned_gradient_x0)

Elemental_pm{Float64}(2, 3, Elemental_em{Float64}[Elemental_em{Float64}(2, [1, 2], [1.0 0.0; 0.0 1.0], Counter_elt_mat(0, 0, 0, 0, 0, 0), false, false, [6.9342354942711e-310, 2.121995815e-314]), Elemental_em{Float64}(2, [2, 3], [1.0 0.0; 0.0 1.0], Counter_elt_mat(0, 0, 0, 0, 0, 0), false, false, [2.1219957915e-314, 6.9340456461372e-310])], sparse(Int64[], Int64[], Float64[], 3, 3), sparse(Int64[], Int64[], Float64[], 3, 3), [[1], [1, 2], [2]], [1, 2, 3])

It can be visualized with the Matrix constructor

Matrix(partitioned_matrix)

3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  2.0  0.0
 0.0  0.0  1.0

The second term of the diagonal accumulates two components of value 1.0 from the two initial element approximations.

The module PartitionedVectors.jl eases the manipulation of the partitioned vectors. In PartitionedStructures.jl, the manipulation of partitioned vectors require the use of routines.

Then you compute the partitioned gradient at x1

partitioned_gradient_x1 = similar(partitioned_gradient_x0)
set_epv!(partitioned_gradient_x1, vector_gradient_element(x1, U))

Elemental_pv{Float64}(2, 3, Elemental_elt_vec{Float64}[Elemental_elt_vec{Float64}([4.0, 1.0], [1, 2], 2), Elemental_elt_vec{Float64}([10.0, 10.0], [2, 3], 2)], [6.9342354972489e-310, 9.3983193583606e-311, 6.9342354660287e-310], [[1], [1, 2], [2]], [1, 2, 3])

and compute the elementwise difference of the partitioned gradients partitioned_gradient_x1 - partitioned_gradient_x0

# copy to avoid side effects on partitioned_gradient_x0
partitioned_gradient_difference = copy(partitioned_gradient_x0)

# apply in place an unary minus to every element gradient
minus_epv!(partitioned_gradient_difference)

# add in place the element vector of partitioned_gradient_x1
# to the corresponding element vector of partitioned_gradient_difference
add_epv!(partitioned_gradient_x1, partitioned_gradient_difference)

# compute the vector y
build_v!(partitioned_gradient_difference)
get_v(partitioned_gradient_difference) == y

true

Then you can define the partitioned quasi-Newton PBFGS update

# apply the partitioned PBFGS update to partitioned_matrix
# and return Matrix(partitioned_matrix)
B_PBFGS = update(partitioned_matrix, partitioned_gradient_difference, s; name=:pbfgs, verbose=true)

3×3 Matrix{Float64}:
 2.75  0.25     0.0
 0.25  4.45909  2.14545
 0.0   2.14545  2.42727

which respects the sparsity structure of ∇²f. Note that there exists a variant of update where s is a partitioned vector similar to partitioned_gradient_difference.

In addition, update() indicates the number of elements: updated, not updated or untouched, as long as the user doesn't set verbose=false. The partitioned update verifies the secant equation

norm(B_PBFGS*s - y) < atol

true

which may also be calculated with

# compute the product partitioned matrix vector
Bs = mul_epm_vector(partitioned_matrix, s)
norm(Bs - y) < atol

true

Other partitioned quasi-Newton approximations

There exist two categories of partitioned quasi-Newton updates. In the first category, each element Hessian $\nabla^2 f_i$ is approximated with a dense matrix, for example: PBFGS. In the second category, each element Hessian $\nabla^2 f_i$ is approximated with a quasi-Newton linear-operator.

Partitioned quasi-Newton operators

Once the partitioned matrix is allocated,

partitioned_matrix_PBFGS = epm_from_epv(partitioned_gradient_x0)
partitioned_matrix_PSR1 = epm_from_epv(partitioned_gradient_x0)
partitioned_matrix_PSE = epm_from_epv(partitioned_gradient_x0)

bools = [true, true, false, true] # false means non convex
partitioned_matrix_PCS = epm_from_epv(partitioned_gradient_x0; convex_vector = bools)

Elemental_pm{Float64}(2, 3, Elemental_em{Float64}[Elemental_em{Float64}(2, [1, 2], [1.0 0.0; 0.0 1.0], Counter_elt_mat(0, 0, 0, 0, 0, 0), true, false, [6.93422568146887e-310, 6.9342352305254e-310]), Elemental_em{Float64}(2, [2, 3], [1.0 0.0; 0.0 1.0], Counter_elt_mat(0, 0, 0, 0, 0, 0), true, false, [0.0, 0.0])], sparse(Int64[], Int64[], Float64[], 3, 3), sparse(Int64[], Int64[], Float64[], 3, 3), [[1], [1, 2], [2]], [1, 2, 3])

you can apply on it any of the four partitioned updates : PBFGS, PSR1, PCS and PSE (by default) :

• PBFGS updates each element approximation with BFGS;
• PSR1 updates each element approximation with SR1;
• PCS updates each element approximation with BFGS if the element is convex and SR1 if it is non convex (designated with bools);
• In PSE, each element approximation is updated with BFGS if it is possible or with SR1 otherwise.

B_PBFGS = update(partitioned_matrix_PBFGS, partitioned_gradient_difference, s; name=:pbfgs)
B_PSR1 = update(partitioned_matrix_PSR1, partitioned_gradient_difference, s; name=:psr1)
B_PCS = update(partitioned_matrix_PSR1, partitioned_gradient_difference, s; name=:pcs)
 # name=:pse by default
B_PSE = update(partitioned_matrix_PSE, partitioned_gradient_difference, s)

3×3 Matrix{Float64}:
 2.75  0.25     0.0
 0.25  4.45909  2.14545
 0.0   2.14545  2.42727

All these approximations satisfy the secant equation as long as every element approximation is updated

norm(mul_epm_vector(partitioned_matrix_PBFGS, s) - y) < atol

true

norm(mul_epm_vector(partitioned_matrix_PSR1, s) - y) < atol

true

norm(mul_epm_vector(partitioned_matrix_PCS, s) - y) < atol

false

norm(mul_epm_vector(partitioned_matrix_PSE, s) - y) < atol

true

Limited-memory partitioned quasi-Newton operators

These operators exist to apply the partitioned quasi-Newton methods to a partially-separable function with large $n_i$, whose element Hessian approximations can't be store as dense matrices. The limited-memory partitioned quasi-Newton operators allocate a quasi-Newton operator LBFGS or LSR1 defined in LinearOperators.jl for each element Hessian approximation. It defines three approximations:

• PLBFGS, each element approximation is a LBFGSOperator;
• PLSR1, each element approximation is a LSR1Operator;
• PLSE, each element approximation may be a LBFGSOperator or a LSR1Operator.

Contrary to the partitioned quasi-Newton operators, each limited-memory partitioned quasi-Newton operators has a different type

partitioned_linear_operator_PLBFGS = eplo_lbfgs_from_epv(partitioned_gradient_x0)
partitioned_linear_operator_PLSR1 = eplo_lsr1_from_epv(partitioned_gradient_x0)
partitioned_linear_operator_PLSE = eplo_lose_from_epv(partitioned_gradient_x0)

Elemental_plo{Float64}(2, 3, Union{Elemental_elo_bfgs{Float64}, Elemental_elo_sr1{Float64}}[Elemental_elo_bfgs{Float64}(2, [1, 2], Linear operator
  nrow: 2
  ncol: 2
  eltype: Float64
  symmetric: true
  hermitian: true
  nprod:   0
  ntprod:  0
  nctprod: 0

, Counter_elt_mat(0, 0, 0, 0, 0, 0), false), Elemental_elo_bfgs{Float64}(2, [2, 3], Linear operator
  nrow: 2
  ncol: 2
  eltype: Float64
  symmetric: true
  hermitian: true
  nprod:   0
  ntprod:  0
  nctprod: 0

, Counter_elt_mat(0, 0, 0, 0, 0, 0), false)], sparse(Int64[], Int64[], Float64[], 3, 3), sparse(Int64[], Int64[], Float64[], 3, 3), [[1], [1, 2], [2]], [1, 2, 3])

The different types simplify the update method, since no argument name is required to determine which update is applied

B_PLBFGS = update(partitioned_linear_operator_PLBFGS, partitioned_gradient_difference, s)
norm(B_PLBFGS * s - y) < atol

true

B_PLSE = update(partitioned_linear_operator_PLSE, partitioned_gradient_difference, s)
norm(B_PLSE * s - y) < atol

true

B_PLSR1 = update(partitioned_linear_operator_PLSR1, partitioned_gradient_difference, s)
norm(B_PLSR1 * s - y) < atol

true

That's it, you have all the tools to implement a partitioned quasi-Newton method. Enjoy!

Tips

The main issue about the definition of partitioned structures is informing the $f_i$ and $U_i$. To address this issue you may want to take a look at ExpressionTreeForge.jl which detects automatically the partially-separable structure from an ADNLPModels or a JuMP model.