skscope: Sparsity-Constrained Optimization

A Layer (from skscope.layer) acts as a decorator for the objective function. Parameters are processed by the layer(s) before being passed to the loss function. This allows you to implement complex constraints or re-parameterizations.

Key features:

• Sequential Composition: You can pass a list of layers to the layers argument in the solve method. Parameters pass through them in order.
• Custom Layers: You can define your own layers by inheriting from skscope.layer.Identity.

Built-in re-parameterization layers:

• NonNegative: Ensures parameters are non-negative.
• LinearConstraint: Adds a linear constraint (e.g., $\sum \theta_i = 1$).
• SimplexConstraint: Constrains parameters to a simplex.
• BoxConstraint: Constrains parameters within a box.
• OffsetSparse: Adds a constant offset to parameters (e.g., making non-selected features equal to $\mu$ instead of zero).

from skscope import ScopeSolver
import skscope.layer as Layer 
from jax import numpy as jnp

X = jnp.array([[1, 2, 3], [4, 5, 6]])
y = jnp.array([1, 2])

def loss(params):
    return jnp.sum((X @ params - y) ** 2)

solver = ScopeSolver(3, 1)

# Example 1: Using OffsetSparse to make non-selected features equal to 1
solver.solve(
    loss,
    layers=[Layer.OffsetSparse(dimensionality=3, offset=1)], 
)

# Example 2: Using multiple layers (LinearConstraint and NonNegative)
solver.solve(
    loss,
    layers=[
        Layer.LinearConstraint(dimensionality=3, coef=jnp.array([[1, 1, 1]])),
        Layer.NonNegative(dimensionality=3),
    ]
)

Eigen Tensor expressions are evaluated (computed) through several mechanisms. Because expressions are lazy, they are not computed until you explicitly trigger evaluation.

1. Assignment to a Tensor type

Assigning an expression to a Tensor, TensorFixedSize, or TensorMap is the most common way to trigger evaluation.

auto t3 = t1 + t2;             // t3 is an Operation (not evaluated)
auto t4 = t3 * 0.2f;           // t4 is an Operation (not evaluated)
Tensor<float, 3> result = t4;  // The operations are evaluated here

If you know the rank and sizes, assigning to a TensorFixedSize is more efficient:

TensorFixedSize<float, Sizes<4, 4, 2>> result = t5;

2. Using the eval() method

Use .eval() to force an intermediate value to be materialized into a temporary tensor. This is useful for:

• Performance: Preventing a complex sub-expression from being re-computed multiple times (e.g., inside a broadcast() or reshape() chain).
• Correctness (Aliasing): Preventing incorrect results when a tensor is used on both sides of an assignment (e.g., Y = Y / Y.sum()...).

Note: The return value of .eval() is itself an Operation, so it does not trigger evaluation immediately. You must still assign it to a Tensor to see results.

// Forces the sum to be computed once before reshaping/broadcasting
Y = Y / (Y.sum(depth_dim).eval().reshape(dims2d).broadcast(bcast));

3. Using TensorRef for partial access

If you only need to access a few specific elements of a large expression, use TensorRef. This avoids materializing the entire result into memory. The expression is evaluated on-the-fly only for the elements you access.

// ref is an Operation, not a Tensor
TensorRef<Tensor<float, 3>> ref = ((t1 + t2) * 0.2f).exp();

// Evaluation happens here for specific elements
float at_0 = ref(0, 0, 0);

// Example of evaluating a chain of operations
auto t3 = t1 + t2;
auto t4 = t3 * 0.2f;
auto t5 = t4.exp();
Tensor<float, 3> result = t5; // Evaluation occurs here

Eigen Tensors support two memory layouts: ColMajor (default) and RowMajor.

Important Constraints:

• Full Support: Currently, only ColMajor is fully supported. Using RowMajor is not recommended.
• Type Specification: Layout is part of the type. If not specified, ColMajor is assumed.
• Compatibility: All operands in a tensor expression must use the same layout. Mixing layouts will cause a compilation error.
• Layout Swapping: You can change the layout using .swap_layout(). Note that swapping the layout also reverses the order of the dimensions. To preserve dimension order while swapping layout, use .swap_layout().shuffle(shuffle) with a permutation array.

// Explicitly declaring ColMajor (default)
Tensor<float, 3, ColMajor> col_major;

// Explicitly declaring RowMajor
TensorMap<Tensor<float, 3, RowMajor>> row_major(data, ...);

// Simple layout swap (reverses dimension order)
// If col_major is (2, 4), result will be (4, 2)
col_major_result = row_major.swap_layout();

// Swap layout while preserving dimension order using shuffle
// Example: swapping (2, 4) to (2, 4) with RowMajor
array<int, 2> shuffle(1, 0);
col_major_result = row_major.swap_layout().shuffle(shuffle);

A Reduction operation returns a tensor with fewer dimensions than the original by applying a reduction operator (like sum() or maximum()) to slices of the input tensor.

Reduction Dimensions

• You specify dimensions to reduce using an array of integers (type <TensorType>::Dimensions).
• Each integer represents the index of a dimension to be reduced.
• The number of elements in the array must not exceed the rank of the input tensor.
• Each dimension index must be unique (duplicates are not removed).
• Performance Tip: Listing dimensions in increasing order may improve execution speed.

Special Case: All Dimensions

If you call a reduction operation without any parameters, it reduces along all dimensions, resulting in a zero-dimension tensor (a scalar).

// Reduce along one dimension (index 1)
Eigen::Tensor<int, 2> a(2, 3);
a.setValues({{1, 2, 3}, {6, 5, 4}});
Eigen::array<int, 1> dims({1});
Eigen::Tensor<int, 1> b = a.maximum(dims);

// Reduce along all dimensions (returns a scalar/0-dim tensor)
Eigen::Tensor<float, 0> scalar = a.sum();