Advanced API

Hessian and solver

The public differentiator, hessian, gradient, matrix, vector, and solver classes expose the numerical protocol beneath the minimizer. This is a generated block-sparse system specialized to YASPS attributes—not a general sparse linear-algebra library.

Generate a Hessian with diff2

The useful way to construct a populated Hessian is second-order differentiation:

from yasps import differentiator

H = differentiator().diff2(
  source=[elasticity, contact],
  target1=targets,
  target2=targets,
  local_targets=[],
  projection_method=2,
  save_intermediate=False,
  separate_hessian_jacobian=False,
  dynamic_instances=False,
)

Every source must be scalar. Passing an empty source list raises ValueError. differentiator.diff1 is declared in the current package but not implemented; use the gradient attached to a computed Hessian for first-order values.

Direct hessian construction

from yasps import hessian

H_empty = hessian(
  wrt=targets,
  local_targets=[],
  dynamic_instances=False,
)

The constructor creates dimensions and empty storage only. It has no energy terms, index kernels, or placement metadata until the differentiator supplies them. Direct construction is therefore useful primarily for generator work, not ordinary assembly.

Dynamic attributes are not valid global Hessian targets.

Add symbolic Hessian terms

H_total = H_elastic + H_contact

Addition concatenates static and dynamic symbolic term metadata. Both operands must have the same target count, order, and attribute hashes. It does not add the current numerical buffers on the host; call compute() on the combined container to assemble its current values.

This is also how the minimizer creates an active Hessian from all nonignored energy requests.

Assemble a Hessian

H.compute()

or, to direct gradient output into an existing compatible layout:

H.compute(local_gradient=g)

The no-return method performs numerical setup and assembly:

1. on first use, generate static sparse indices and coordinate compression;
2. initialize dynamic sparse metadata when dynamic terms exist;
3. on later dynamic assemblies, refresh changing coordinates;
4. launch placement/reorder and numerical kernels;
5. accumulate Hessian block values, the scalar diagonal, dense diagonal blocks, and gradient values.

Calling H.compute() directly does not construct the inverse diagonal blocks required by the stock PCG preconditioner. minimizer.computeNumericValue adds that step. Use the minimizer to prepare a Hessian before a direct solver.computeSolution call.

Sparse-index lifecycle

The public methods below exist mainly for the generator protocol:

Ordinary applications should call compute() and let it select the correct path. Calling these manually out of order can desynchronize positions, reorder maps, block counts, and numerical buffers.

Block-sparse representation

YASPS groups blocks by their (rows, cols) dimensions. Static and dynamic energy terms use parallel storage:

The matrix base also exposes camelCase aliases such as blockDimensions and blocksFlattened for compatibility.

Diagonal and preconditioner storage

The minimizer builds an inversion kernel specialized to the distinct target element sizes, target instance counts, and block offsets. It reuses that kernel while the target layout stays stable and recomputes the inverse values after each numerical assembly.

Generator metadata

Properties including global_hessians, global_jacobians, indices_kernels, block_indices_gpu, placement_reorder_kernels, and their dynamic forms are public because the Cython subsystems exchange them. Their shapes are part of the code-generation implementation rather than a stable application data model.

Construct and inspect a gradient

from yasps import gradient

g = gradient(
  wrt=targets,
  hessian=H,
)
g.compute()

g.compute() takes no parameters and returns None. With a parent Hessian it calls H.compute(g); without one it leaves the initialized zero vector unchanged.

Segment sizes are numInstances × rows × cols. Segments alias value; later assemblies refill the same storage.

Prepare a direct solver call

The safest source of a direct-solve Hessian is a minimizer, because it assembles both numerical values and inverse diagonal blocks:

import pycuda.gpuarray as gpuarray
from yasps import minimizer, solver

engine = minimizer()
engine.addEnergy(
  energy,
  targets=targets,
  projection_method=2,
)
engine.addWrt(targets)

H = engine.computeNumericValue()
initial_guess = gpuarray.zeros_like(engine.gradient)

If H is None, every energy is ignored and there is no active system to solve.

Call the PCG solver

pcg = solver()
status = pcg.computeSolution(
  active_hessian=H,
  wrt=engine.wrt,
  gradient_object=H.gradient,
  initial_guess=initial_guess,
  tolerance=1e-3,
  maxIterations=20000,
)

flat_direction = pcg.solution

The solver owns reusable work vectors and one generated solverKernel association. On each call it clears the work vectors and its solution, updates the block-dimension list, then invokes the kernel.

status < 0 means the solve did not converge. Exact negative codes are not a documented API. pcg.solution still contains the best result produced.

The solution is a single flattened buffer. Use H.gradient.gradient_segments_start_cpu, or the known target segment sizes, to construct target-aligned views.

Reset solver state

pcg.reset()

This no-argument method returns None. It discards the compiled solver-kernel association, all PCG work buffers, and the solution buffer. Repeated solves with a stable layout should not call reset; they reuse those allocations.

Reset after changing to an incompatible target or block-dimension layout. minimizer.addWrt resets its owned solver automatically.

Empty-system behavior

If the gradient has zero elements, computeSolution returns 0. If no active Hessian is supplied, it returns 0 and leaves or clears the solution as appropriate. A normal nonempty solve requires a non-None Hessian.

The vector class

from yasps import vector

v = vector(size=number_of_scalars)

The current CUDA implementation stores float64 values.

v.updateValue(numpy_array_or_gpuarray_or_vector)
w = v + other
w = v - other
w = -v
w = 2.0 * v

updateValue(new_value) requires exactly size elements. NumPy and non-float64 GPU inputs are converted to float64. Binary operands must have equal sizes, and arithmetic creates a new result vector.

resize(new_size) updates the recorded size only; it does not reallocate the underlying GPU buffer. Treat it as low-level metadata, not a safe general resize operation.

The base matrix class

from yasps import matrix

M = matrix(
  rows=0,
  cols=0,
)

matrix is a storage base. Its matVecProduct and matVecProductInPlace methods currently validate dimensions but do not implement block-sparse multiplication. Consequently, matrix * vector is not a usable numerical Hessian-vector product.

PCG multiplication is generated by solverKernel and operates directly on the Hessian’s static and dynamic block buffers.

Safe inspection versus mutation

Reading numerical blocks, gradients, diagonals, and segment offsets is useful and intentional. Mutating block positions, counts, reorder metadata, or generator kernel lists requires updating every related component:

sparse coordinates
  → compression and placement maps
  → numerical block ordering
  → diagonal-block preconditioner
  → solver block-dimension protocol

For custom numerical experiments, preserve those structures and replace only the solve that consumes an assembled minimizer Hessian.