High Order Prior: private API

PortfolioOptimisers.block_vec_pq — Function

block_vec_pq(A::MatNum, p::Integer, q::Integer)

Block vectorisation operator.

block_vec_pq transforms a matrix A into a block vectorised form, partitioning A into blocks of size (p, q) and writing each vectorised block as one row. This is useful for higher-order moment computations and tensor manipulations in portfolio analytics.

Mathematical definition

Partition $\mathbf{A}$ into $m \times n$ blocks $\mathbf{A}_{ij}$ of size $p \times q$. The block vectorisation writes each block as one row, in block-column order:

\[\mathcal{V}_{p,q}(\mathbf{A}) = \begin{bmatrix} \mathrm{vec}(\mathbf{A}_{11})^\intercal \\ \vdots \\ \mathrm{vec}(\mathbf{A}_{m1})^\intercal \\ \mathrm{vec}(\mathbf{A}_{12})^\intercal \\ \vdots \\ \mathrm{vec}(\mathbf{A}_{mn})^\intercal \end{bmatrix}\,.\]

Where:

$\mathbf{A}_{ij}$: block $(i, j)$ of $\mathbf{A}$, holding rows $(i-1)p+1$ to $ip$ and columns $(j-1)q+1$ to $jq$.
$\mathrm{vec}$: column-major vectorisation.
$m = \mathrm{size}(\mathbf{A}, 1) / p$, $n = \mathrm{size}(\mathbf{A}, 2) / q$: the block counts.

Row $(j-1)m + i$ of the result is $\mathrm{vec}(\mathbf{A}_{ij})^\intercal$, so the block column index runs slowest. A square partition, $m = n$ and $p = q$, cannot separate this order from the one that runs the block row index slowest, because the two differ by a permutation that is the identity there.

Algorithm

Read size(A) into mp and nq, and check both divisibility conditions.
Take the block counts m = mp ÷ p and n = nq ÷ q.
Allocate A_vec, of size (m * n, p * q).
For each block column j, build Aj, whose i-th row is the vectorisation of block (i, j) of A.
Write Aj into rows j * m + 1 to (j + 1) * m of A_vec.

Arguments

A: Input matrix of size (m * p, n * q), where m and n are integers.
p: Number of rows in each block.
q: Number of columns in each block.

Validation

size(A, 1) must be an integer multiple of p.
size(A, 2) must be an integer multiple of q.

Returns

A_vec::Matrix: Block vectorised matrix of size (m * n, p * q).

Examples

julia> A = [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16];julia> PortfolioOptimisers.block_vec_pq(A, 2, 2)4×4 Matrix{Int64}:  1   5   2   6  9  13  10  14  3   7   4   8 11  15  12  16

Related

dup_elim_sum_matrices

References

[5] D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025). Appendix A.1, Equation A.14.

source

PortfolioOptimisers.elimination_matrix — Function

elimination_matrix(n::Int, diag::Bool = true)

Construct the elimination matrix for a symmetric matrix of size n × n.

The elimination matrix L extracts the unique (lower triangular) elements of a symmetric matrix. Used internally in coskewness and cokurtosis computation.

Mathematical definition

$\mathbf{L}_n$ is defined by the identity it applies, for every $\mathbf{A}$ of size $n \times n$:

\[\mathbf{L}_n \mathrm{vec}(\mathbf{A}) = \mathrm{vech}(\mathbf{A})\,.\]

With diag = false it drops the diagonal as well:

\[\mathbf{L}_n^{-} \mathrm{vec}(\mathbf{A}) = \mathrm{vech}^{-}(\mathbf{A})\,.\]

Where:

$\mathrm{vec}$: column-major vectorisation, of length $n^2$.
$\mathrm{vech}$: half-vectorisation, the lower triangle read column by column, of length $n(n+1)/2$.
$\mathrm{vech}^{-}$: the strictly lower triangle read the same way, of length $n(n-1)/2$.

The identity holds for any square $\mathbf{A}$, symmetric or not, because $\mathbf{L}_n$ only reads the lower triangle. Dropping the diagonal removes the $n$ rows that read it, so the row count falls from $n(n+1)/2$ to $n(n-1)/2$ while the column count stays $n^2$.

Algorithm

Take nsq = n^2.
Read diag, and set the row count m, the column range rg and the offset b from it: m = n(n+1)/2, rg = 1:n and b = 0 when diag is true, and m = n(n-1)/2, rg = 2:n and b = 1 otherwise.
Fill v, whose r-th entry is the position in $\mathrm{vec}(\mathbf{A})$ of the r-th entry of the half-vectorisation. b carries the offset that skips the entries above the diagonal, and — under diag = false — the diagonal entry too.
Return the sparse matrix carrying a one at each (r, v[r]), of size m × nsq.

Arguments

n: Size of the symmetric matrix.
diag: Whether to include the diagonal elements.

Returns

Sparse elimination matrix.

Related

duplication_matrix
summation_matrix
dup_elim_sum_matrices: builds this matrix and its two siblings in one walk.

References

[5] D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025). Appendix A.2, Equation A.26.

source

PortfolioOptimisers.summation_matrix — Function

summation_matrix(n::Int, diag::Bool = true)

Construct the summation matrix for a symmetric matrix of size n × n.

The summation matrix S adds up contributions from both triangular halves of a symmetric matrix. Used internally in coskewness and cokurtosis computation.

Mathematical definition

$\mathbf{S}_n$ is the elimination matrix reweighted by the multiplicity that $\mathbf{D}_n$ restores:

\[\mathbf{S}_n = \mathbf{D}_n^\intercal \mathbf{D}_n \mathbf{L}_n\,,\]

so that, for every $\mathbf{A}$ of size $n \times n$,

\[\boldsymbol{1}^\intercal \mathbf{S}_n \mathrm{vec}(\mathbf{A}) = \boldsymbol{1}^\intercal \mathrm{vec}(\mathbf{A})\,.\]

Where:

$\mathbf{D}_n$: the duplication matrix of duplication_matrix.
$\mathbf{L}_n$: the elimination matrix of elimination_matrix.
$\boldsymbol{1}$: vector of ones of the length its neighbour needs.

$\mathbf{S}_n$ reads the lower triangle and weights each entry by the number of places it occupies in $\mathrm{vec}(\mathbf{A})$: one for a diagonal entry and two for an off-diagonal one. The sum identity above follows, and it holds for any square $\mathbf{A}$, symmetric or not.

The same construction with diag = false gives $\mathbf{S}_n^{-} = (\mathbf{D}_n^{-})^\intercal \mathbf{D}_n^{-} \mathbf{L}_n^{-}$, which weights every one of its $n(n-1)/2$ rows by two. Its sum identity therefore reaches the off-diagonal entries alone, $\boldsymbol{1}^\intercal \mathbf{S}_n^{-} \mathrm{vec}(\mathbf{A}) = \boldsymbol{1}^\intercal \mathrm{vec}(\mathbf{A} - \mathrm{diag}(\mathbf{A}))$.

Algorithm

The body builds the product of the definition directly, without forming $\mathbf{D}_n$.

Take nsq = n^2. Read diag, and set the row count m, the column range rg and the offset b from it, exactly as elimination_matrix does.
Walk the columns in rg. Write into v1 the $\mathrm{vec}$ position of every entry of the half-vectorisation, and into v2 and rows2 the $\mathrm{vec}$ position and the half-vectorisation row of each strictly lower entry.
Drop the zero entries of v1, v2 and rows2, which are the slots the walk never filled.
When diag is true, return the sum of two sparse matrices: one carrying a one at each (r, v1[r]), which is $\mathbf{L}_n$, and one carrying a one at each (rows2[k], v2[k]), which adds the second unit to every off-diagonal row.
When diag is false, every row is off-diagonal, so return the sparse matrix carrying a two at each (r, v1[r]). v2 and rows2 go unread on this branch.

Arguments

n: Size of the symmetric matrix.
diag: Whether to include the diagonal elements.

Returns

Sparse summation matrix.

Related

duplication_matrix
elimination_matrix
dup_elim_sum_matrices: builds this matrix and its two siblings in one walk, through the product of the definition rather than by the construction above.

References

[5] D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025). Appendix A.2, Equation A.27.
[30] D. Cajas. Convex Optimization of Portfolio Kurtosis. Available at SSRN 4202967 (2022).

source

PortfolioOptimisers.dup_elim_sum_matrices — Function

dup_elim_sum_matrices(n::Int)

Construct duplication, elimination, and summation matrices for symmetric matrix vectorisation.

dup_elim_sum_matrices returns the duplication matrix D, elimination matrix L, and summation matrix S for symmetric matrices of size N × N. These matrices are used in higher-order moment computations, tensor manipulations, and efficient vectorisation of symmetric matrices in portfolio analytics.

The three are the diag = true matrices of duplication_matrix, elimination_matrix and summation_matrix, built in one walk of the columns rather than in three.

Mathematical definition

For every $\mathbf{A}$ of size $n \times n$:

\[\begin{align} \mathbf{D}_n \mathrm{vech}(\mathbf{A}) &= \mathrm{vec}(\mathbf{A})\,, \\ \mathbf{L}_n \mathrm{vec}(\mathbf{A}) &= \mathrm{vech}(\mathbf{A})\,, \\ \mathbf{S}_n &= \mathbf{D}_n^\intercal \mathbf{D}_n \mathbf{L}_n\,. \end{align}\]

Where:

$\mathrm{vec}$: column-major vectorisation, of length $n^2$.
$\mathrm{vech}$: half-vectorisation, the lower triangle read column by column, of length $n(n+1)/2$.

The third line makes $\boldsymbol{1}^\intercal \mathbf{S}_n \mathrm{vec}(\mathbf{A}) = \boldsymbol{1}^\intercal \mathrm{vec}(\mathbf{A})$: $\mathbf{S}_n$ reads the lower triangle and weights each off-diagonal entry by the two places it occupies in $\mathrm{vec}(\mathbf{A})$.

Algorithm

Check that n is positive.
Take m = n(n+1)/2 and nsq = n^2.
Walk the columns once, filling v1 and v2 together. v1 is the column index vector of duplication_matrix and v2 is that of elimination_matrix.
Build d, the sparse matrix carrying a one at each (r, v1[r]), of size nsq × m.
Build l, the sparse matrix carrying a one at each (r, v2[r]), of size m × nsq.
Build s as the product transpose(d) * d * l of the definition above, rather than by the direct construction of summation_matrix.

Arguments

n: Size of the symmetric matrix (integer).

Validation

n > 0.

Returns

D::SparseMatrixCSC{Int64, Int64}: Duplication matrix (n^2 × m), where m = n(n+1)/2.
L::SparseMatrixCSC{Int64, Int64}: Elimination matrix (m × n^2).
S::SparseMatrixCSC{Int64, Int64}: Summation matrix (m × n^2).

Examples

julia> D, L, S = PortfolioOptimisers.dup_elim_sum_matrices(3);julia> D9×6 SparseArrays.SparseMatrixCSC{Int64, Int64} with 9 stored entries: 1  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  1  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  1  ⋅  ⋅  ⋅ ⋅  1  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  1  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  1  ⋅ ⋅  ⋅  1  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  1  ⋅ ⋅  ⋅  ⋅  ⋅  ⋅  1julia> L6×9 SparseArrays.SparseMatrixCSC{Int64, Int64} with 6 stored entries: 1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1julia> S6×9 SparseArrays.SparseMatrixCSC{Int64, Int64} with 6 stored entries: 1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  2  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  2  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  ⋅  2  ⋅  ⋅  ⋅ ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1

Related

block_vec_pq
duplication_matrix
elimination_matrix
summation_matrix

References

[5] D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025). Appendix A.2, Equations A.25 to A.27.

source

PortfolioOptimisers.duplication_matrix — Function

duplication_matrix(n::Int, diag::Bool = true)

Construct the duplication matrix for a symmetric matrix of size n × n.

The duplication matrix D maps the vech (half-vectorisation) of a symmetric matrix to its full vec. Used internally in coskewness and cokurtosis computation.

Mathematical definition

$\mathbf{D}_n$ is defined by the identity it restores, for every symmetric $\mathbf{A}$ of size $n \times n$:

\[\mathbf{D}_n \mathrm{vech}(\mathbf{A}) = \mathrm{vec}(\mathbf{A})\,.\]

With diag = false the half-vectorisation drops the diagonal, and the identity restores the hollow matrix:

\[\mathbf{D}_n^{-} \mathrm{vech}^{-}(\mathbf{A}) = \mathrm{vec}(\mathbf{A} - \mathrm{diag}(\mathbf{A}))\,.\]

Where:

$\mathrm{vec}$: column-major vectorisation, of length $n^2$.
$\mathrm{vech}$: half-vectorisation, the lower triangle read column by column, of length $n(n+1)/2$.
$\mathrm{vech}^{-}$: the strictly lower triangle read the same way, of length $n(n-1)/2$.
$\mathrm{diag}(\mathbf{A})$: the diagonal of $\mathbf{A}$ held as a matrix.

Each row of $\mathbf{D}_n$ carries exactly one entry, so the matrix selects rather than sums. A diagonal entry of $\mathbf{A}$ is selected once and an off-diagonal entry twice, which is what makes $\mathbf{D}_n^\intercal \mathbf{D}_n$ the weight matrix that summation_matrix applies.

Algorithm

Take m = n(n+1)/2 and nsq = n^2.
Fill v, whose r-th entry is the position in $\mathrm{vech}(\mathbf{A})$ of the entry that row r of $\mathrm{vec}(\mathbf{A})$ holds. The inner loops walk the strictly upper part of a column first, then its lower part.
When diag is true, return the sparse matrix carrying a one at each (r, v[r]), of size nsq × m.
When diag is false, count how often each position occurs in v, giving counts.
Keep the positions that occur more than once — the off-diagonal ones — and renumber them from one, giving cols.
Keep every row of $\mathrm{vec}(\mathbf{A})$ whose position survives step 5, giving filtered_rows and filtered_cols.
Return the sparse matrix carrying a one at each kept pair, of size nsq × n(n-1)/2.

Arguments

n: Size of the symmetric matrix.
diag: Whether to include the diagonal elements.

Returns

Sparse duplication matrix.

Related

elimination_matrix
summation_matrix
dup_elim_sum_matrices: builds this matrix and its two siblings in one walk.

References

[5] D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025). Appendix A.2, Equation A.25.

source

PortfolioOptimisers.dup_elim_sum_view — Method

dup_elim_sum_view(args...)

Answer with three nothings when the first argument is not a matrix.

The fallback of dup_elim_sum_view. It builds nothing and reads none of its arguments; the matrix method is the one that calls dup_elim_sum_matrices.

Its two call sites are in port_opt_view on a HighOrderPrior, which passes pr.kt as the first argument. No estimator in the library builds a carrier that reaches this method: sk and V travel together, and kt, L2 and S2 do too, so a carrier holding any of them holds kt. A hand-built carrier can — D2 is the one moment field the constructor accepts on its own — and it is the case this method answers.

This is a varargs fallback, so it also answers any call whose argument count the matrix method does not take. dup_elim_sum_view(M, n) reaches the matrix method for a matrix M; dup_elim_sum_view(M, n, extra) reaches this one.

Algorithm

Return (nothing, nothing, nothing), reading no argument.

Arguments

args...: Any arguments. None is read.

Returns

(nothing, nothing, nothing)::Tuple{Nothing, Nothing, Nothing}.

Related

duplication_matrix
elimination_matrix
summation_matrix

source

PortfolioOptimisers.dup_elim_sum_view — Method

dup_elim_sum_view(
    _::AbstractMatrix{<:Union{var"#s136", var"#s53"} where {var"#s136"<:Number, var"#s53"<:AbstractJuMPScalar}},
    n
) -> Tuple{SparseArrays.SparseMatrixCSC{Int64, Int64}, SparseArrays.SparseMatrixCSC{Int64, Int64}, SparseArrays.SparseMatrixCSC{Int64, Int64}}

Compute duplication, elimination, and summation matrices at the dimension the caller names.

Overload of dup_elim_sum_view for a matrix first argument. The matrix is read for dispatch alone, and the dimension is the second argument n, not size of the matrix. port_opt_view on a HighOrderPrior relies on that: it passes the carrier's full N^2 × N^2 cokurtosis and the asset count of the subproblem, so the three matrices come back rebuilt at the smaller dimension rather than cut from the larger ones.

Algorithm

Forward n to dup_elim_sum_matrices, and return its three matrices.

Arguments

The first argument: any matrix. It selects this method and is not read.
n: Size of the symmetric matrix the three matrices are built for.

Returns

(D, L, S)::Tuple{SparseMatrixCSC, SparseMatrixCSC, SparseMatrixCSC}: The three matrices of dup_elim_sum_matrices at dimension n.

Related

duplication_matrix
elimination_matrix
summation_matrix

source

PortfolioOptimisers.assemble_high_order_prior — Function

assemble_high_order_prior(
    pe::HighOrderPriorEstimator,
    pr::AbstractPriorResult,
    kt::Union{Nothing, AbstractMatrix{<:Union{var"#s136", var"#s53"} where {var"#s136"<:Number, var"#s53"<:AbstractJuMPScalar}}},
    sk::Union{Nothing, AbstractMatrix{<:Union{var"#s136", var"#s53"} where {var"#s136"<:Number, var"#s53"<:AbstractJuMPScalar}}},
    V
) -> HighOrderPrior{var"#s185", _A, _B, _C, _D, _E, _F, _G, Nothing} where {var"#s185"<:AbstractPriorResult, _A, _B, _C, _D, _E, _F, _G}

Assembles a HighOrderPrior from a low order result and the co-moments fitted beside it.

The tail of every HighOrderPriorEstimator fit, written once: the batch method reaches it after running the two co-moment verbs over the caller's matrix, and the read-out of a folded estimator reaches it after reading the same two off their states. The duplication-elimination matrices a consumer needs depend on which co-moments are present and on nothing else, so the rule lives here rather than in each caller.

D2, L2 and S2 are sized from the full width of pr.X, which is what the expanded tensors carry.

Algorithm

Build D2, L2 and S2 where both co-moments are present, and L2 and S2 alone where only the cokurtosis is.
Assemble the HighOrderPrior, carrying pe.ske.mp where a coskewness tensor was fitted.
Refuse a carrier whose blocks do not agree on the Coverage Universe, with assert_matched_coverage.

Arguments

pe: High order prior estimator.
pr: The low order prior result the embedded estimator answered.
kt: The square cokurtosis matrix, or nothing.
sk: The coskewness tensor, or nothing.
V: The coskewness view matrix, or nothing.

Validation

Every block of the carrier agrees on the Coverage Universe. An error is thrown otherwise.

Returns

hop::HighOrderPrior: The assembled result.

Related

source

PortfolioOptimisers.comoment_investable — Function

comoment_investable(
    sk::Union{Nothing, AbstractMatrix{<:Union{var"#s136", var"#s53"} where {var"#s136"<:Number, var"#s53"<:AbstractJuMPScalar}}},
    kt::Union{Nothing, AbstractMatrix{<:Union{var"#s136", var"#s53"} where {var"#s136"<:Number, var"#s53"<:AbstractJuMPScalar}}},
    N::Integer
) -> BitVector

Reads the per-asset diagonal of every co-moment a HighOrderPrior holds, and answers where all of them are finite.

The higher-order half of the Investable Mask. A coskewness tensor's per-asset diagonal is sk[i, (i - 1) * N + i], which is asset i's third central moment, and a cokurtosis matrix's is kt[j, j] at the pair column j = (i - 1) * N + i, which is its fourth. Both are finite exactly where the fit answered for that asset, which is the same rule investable_mask reads off the diagonal of sigma. A moment the carrier does not hold constrains nothing, so a nothing reads as every asset admitted.

Arguments

sk: The coskewness tensor, assets × assets², or nothing.
kt: The square cokurtosis matrix, assets² × assets², or nothing.
N: Number of assets.

Returns

msk::BitVector: true at every asset whose higher-order moments are finite.

Related

source

PortfolioOptimisers.assert_matched_coverage — Function

assert_matched_coverage(pr::HighOrderPrior)

Says so, once, when the higher orders of a fitted HighOrderPrior cover fewer assets than its low order block does.

A CoveragePolicy is a field of one estimator, so a caller may set one on pe.pe's mean and covariance and leave pe.ske and pe.kte on the Coverage Universe. That configuration is legal and well defined — the low orders answer an asset that lists inside the window and the higher orders do not — but the Investable Mask then narrows back to the Coverage Universe, and the caller has bought nothing where the panel is gappiest. It is only the silence that is refused, exactly as scenario_fill refuses it for a zero-filled scenario.

The check runs once, at the fit, rather than in investable_mask, which a fold loop calls at every optimiser entry.

Arguments

pr: The fitted high order prior.

Validation

Every asset the low order block holds is held by the higher orders too. A warning naming the assets is emitted otherwise.

Returns

nothing.

Related

source

References

[5]: D. Cajas. Advanced Portfolio Optimization: A Cutting-edge Quantitative Approach (Springer Nature Switzerland, 2025).
[30]: D. Cajas. Convex Optimization of Portfolio Kurtosis. Available at SSRN 4202967 (2022).