A PartitionedVector <: DenseVector <: AbstractVector wraps a PartitionedStructures.Elemental_pv, to make JuliaSmoothOptimizers modules able to exploit the partially separable structure. For now, PartitionedVector is the keystone to define PartiallySeparableNLPModels, and then to consequently:
- replace
Vector in JSOSolvers.jl; - replace
Vector in KrylovSolvers.jl; - make
LinearOperators relying only on PartitionedVectors viable; - fit the NLPModels interface (ex:
NLPModels.obj(nlp::PartiallySeparableNLPModel, x::PartitionedVector) see PartiallySeparableNLPModels.jl for more details).
Partitioned vectors usages are related to partially separable problems:
\[ 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,\]
whose inherit partitioned derivatives:
\[ \nabla f(x) = \sum_{i=1}^N U_i^\top \nabla f_i (U_i x) \quad \nabla^2 f(x) = \sum_{i=1}^N U_i^\top \nabla^2 f_i (U_i x) U_i,\]
where the set of $U_i$ gathers all the information related to the partitioned structure. We consider $U_i$ a linear application selecting the variable that parametrize the $i$-th element function $f_i$, each $U_i$ is defined from a vector of integers. The partially separable function
\[f(x) = f_1(x_1,x_2,x_3,x_4) + f_2(x_3,x_4,x_5,x_6) + f_3(x_5,x_6,x_7) + f_4(x_5,x_6,x_8)\]
get the following $U_i$:
U1 = [1, 2, 3, 4]
+U2 = [3, 4, 5, 6]
+U3 = [5, 6, 7]
+U4 = [5, 6, 8]
+Uis = [U1, U2, U3, U4]
4-element Vector{Vector{Int64}}:
+ [1, 2, 3, 4]
+ [3, 4, 5, 6]
+ [5, 6, 7]
+ [5, 6, 8]
However, before defining a PartitionedVector, you must understand the two usages of PartitionedVectors:
- usage 1: store each element contribution $\hat{v}_i$ independently and aggregate them to form $v = \sum_{i=1}^N U_i^\top \hat{v}_i$
In optimization methods, you use it to store $\nabla f_i (U_i x)$, $y_i = \nabla f_i (U_i x_{k+1}) - \nabla f_i (U_i x_k)$ or the result of a partitioned hessian-vector product which is a partitioned vector:
\[ \nabla^2 f(x_k) s \approx B_k s = (\sum_{i=1}^N U_i^\top B_{k,i} U_i ) s = \sum_{i=1}^N U_i^\top (B_{k,i} U_i s) .\]
- usage 2: represent simultaneously a vector $x \in \mathbb{R}^n$ and the application of every $U_i$ onto $x$: $U_i x, \forall i$.
By construction, the elements parametrized by the same variables (for exemple U1 and U3 are parametrized by the third variable) share the same values. In optimization methods, it allows to store the current point $x_k$ or step $s_k$, which always comes in handy to evaluate $f_i(U_i x), \nabla f_i(U_i x)$ or $B_{k,i} U_i s$.
Any methods exploiting partially separable concepts will have to manipulate both usages at the same time, in particular the solvers from JSOSolvers.jl and Krylov.jl. By default, you can define a PartitionedVector just from Uis as
using PartitionedVectors
+pv = PartitionedVector(Uis)
PartitionedVector{Float64}
+[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
which will be of usage 1. If you want a vector for usage 2 you have to set simulate_vector optional argument to true
pv_vec = PartitionedVector(Uis; simulate_vector=true)
PartitionedVector{Float64}
+[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
The distinction between both is about how they form a Vector from them. For example:
pv .= pv_vec
+pv == pv_vec
which means that elements from pv are of same value than those of pv_vec. However, Vector(pv) relies on build! satisfying the specified usage and consequently
Vector(pv) != Vector(pv_vec)
You can set a PartitionedVector of usage 2 from a Vector with
set!(pv_vec, rand(length(pv_vec)))
PartitionedVector{Float64}
+[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Warning: the application set!() on PartitionedVector of usage 1 doesn't have sense, and produces an error.
PartitionedVectors.jl specify several methods from various modules. For operations that are not in place, the result will take usage 1.
Warning: you have to be careful when you mix both usages in a single operation, because the result may have no meaning. Keep in mind the result you expect for Vector(pv).
Base.jl:
- elementary operations
+, -, *, == for PartitionedVectors.
pv + pv == 2 * pv
+pv - pv == 0 * pv
true
Moreover, it supports the broadcast for the same operations (unfortunately not in place for now).
pv .+ pv == 2 .* pv
+pv .- pv == 0 .* pv
true
Special case when both usages get the same Vector
pv .= 0
+pv_vec .= pv
+Vector(pv) == Vector(pv_vec)
true
- supports methods related to
AbstractArrays such as: get_index, set_index! refers to the i-th element vector of pv
pv .= 1
+pv[1]
pv2 = pv - 0.5 * pv
+pv[1] = pv2[1]
+pv[1] == pv[2]
firstindex, lastindex, size, length:
firstindex(pv) == 1
+lastindex(pv) == 4
+size(pv) == (4,)
+length(pv) == 8
note that size and length differ, length refers to the biggest index (e.g. $n$) while size refers to the number of elements $N$.pvcopy = copy(pv)
+pvsimilar = similar(pv)
LinearAlgebra: both dot and norm rely on build!(::PartitionedVector) before applying dot or norm to the resulting Vector
using LinearAlgebra
+dot(pv,pv) ≈ norm(pv)^2
true
In addition, PartitionedVectors.jl defines a:
LinearOperator from a PartitionedMatrix (see PartitionedStructures.jl) relying on PartitionedVectors
using PartitionedStructures, LinearOperators
+epm = PartitionedStructures.identity_epm(Uis)
+lo = LinearOperators.LinearOperator(epm)
+lo * pv
PartitionedVector{Float64}
+[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Note: Matrix(lo) will produce an error, since the default implementation assumes a complete Vector-like behaviour.
- dedicated
CGSolver from Krylov.jl to solve a partitioned linear system (from a partitioned LinearOperator).
using Krylov # 0.9.0
+solver = Krylov.CgSolver(pv)
+
+pv_gradient = similar(pv)
+(N,) = size(pv)
+for i = 1:N
+ nie = pv_gradient[i].nie #size of a element vector
+ pv_gradient[i] = rand(nie)
+end
+
+Krylov.solve!(solver, lo, -pv_gradient)
┌──────────────────┬──────────────────────────┬──────────────────┐
+│ CgSolver│ nrows: 8│ ncols: 8│
+├──────────────────┼──────────────────────────┼──────────────────┤
+│Precision: Float64│ Architecture: GPU│Storage: 163 bytes│
+├──────────────────┼──────────────────────────┼──────────────────┤
+│ Attribute│ Type│ Size│
+├──────────────────┼──────────────────────────┼──────────────────┤
+│ m│ Int64│ 8 bytes│
+│ n│ Int64│ 8 bytes│
+│ Δx│PartitionedVector{Float64}│ 16 bytes│
+│ x│PartitionedVector{Float64}│ 16 bytes│
+│ r│PartitionedVector{Float64}│ 16 bytes│
+│ p│PartitionedVector{Float64}│ 16 bytes│
+│ Ap│PartitionedVector{Float64}│ 16 bytes│
+│ z│PartitionedVector{Float64}│ 16 bytes│
+│ warm_start│ Bool│ 1 byte│
+│ stats│Krylov.SimpleStats{Float64}│ 50 bytes│
+└──────────────────┴──────────────────────────┴──────────────────┘
+
+Simple stats
+ niter: 3
+ solved: true
+ inconsistent: false
+ residuals: []
+ Aresiduals: []
+ κ₂(A): []
+ status: solution good enough given atol and rtol
+