Add new initializer: hyper_init

[sm2909]

Fixes #303

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the contributor guidelines, in particular the SciML Style Guide and COLPRAC.
• Any new documentation only uses public API

Implementation Details

The hyper_init function executes the following steps:

1. Sampling: Samples $N$ points on a Poincaré ball.
2. Distance Calculation: Calculates hyperbolic distances between all sampled points, populating an $N \times N$ kernel.
3. Weighting: Applies $\sigma$ on the kernel to establish geometry-aware kernel weights.
4. Sparsification: Sparsifies the kernel by selecting the top $k$ values in each row.
5. Scaling: Applies spectral scaling to finalize the matrix.

Refactoring & Helper Functions

To keep hyper_init clean, the logic is modularized into three helper functions (taking inspiration from the original Proof of Concept):

• hyperbolic_distance_poincare
• sample_poincare_points
• topk_row_sparsify

Mathematical Assumptions

A simplifying assumption is used when sampling the radius $\rho$ from the Poincaré ball. To match the approximation used in the PoC, the radius is sampled linearly:

rho = rand(rng) * R_max

Instead of the exact inverse transform sampling method:

u = rand(rng)
rho = acosh(u * (cosh(R_max) - 1) + 1)

Questions:

1. Method Signature: Could you please take a look at the method signature? If the kwargs need more clarity or better naming conventions, I am happy to adjust them.
2. Verification & Testing: How can I formally verify that the initialization values produced are mathematically accurate for a hyperbolic init? The procedural implementation matches the PoC from my understanding. If there is a method to test this, I am happy to test it.

[MartinuzziFrancesco]

remove #

[MartinuzziFrancesco]

better to avoid single letter kwargs

[MartinuzziFrancesco]

Thanks for looking into this! The implementation seems to be ok, although I will need to double check against the python code. There are a couple of minor changes to tackle in the meantime, I left some comments

[MartinuzziFrancesco]

out of curiosity, have you tried to compare it with the python implementation?

[sm2909]

Yes, I have written the current implementation taking direct inspiration from the python code.

The only major difference is:
The python code does sample_poincare_points only on a 2 dimensional ball, whereas the implementation in this PR samples from a general poincare_dim poincare ball.

[MartinuzziFrancesco]

could you try and post the results for a couple of matrices with the python and julia versions?

[sm2909]

Python version:
1.

>>> model = HypER3D(
...     reservoir_size=5, sigma=1.0, top_k=2,
...     spectral_radius=1.0, input_scale=0.2, leaking_rate=0.8, seed=42
... )
... model.build()
... 
>>> model.W
array([[0.        , 0.        , 0.0234523 , 0.17650021, 0.        ],
       [0.        , 0.        , 0.        , 0.25385294, 0.03823212],
       [0.        , 0.        , 0.        , 0.33375124, 0.05703526],
       [0.        , 0.        , 0.33375124, 0.        , 0.9220911 ],
       [0.        , 0.        , 0.05703526, 0.9220911 , 0.        ]])

>>> model = HypER3D(
...     reservoir_size=5, sigma=1.0, top_k=2,
...     spectral_radius=1.0, input_scale=0.2, leaking_rate=0.8, seed=41
... )
... model.build()
... 
>>> model.W
array([[0.        , 0.        , 0.        , 0.48833225, 0.16386675],
       [0.        , 0.        , 0.        , 0.63223134, 0.12939502],
       [0.08204594, 0.        , 0.        , 0.03723776, 0.        ],
       [0.        , 0.63223134, 0.        , 0.        , 0.68339259],
       [0.16386675, 0.        , 0.        , 0.68339259, 0.        ]])

>>> model = HypER3D(
...     reservoir_size=5, sigma=1.0, top_k=2,
...     spectral_radius=1.0, input_scale=0.2, leaking_rate=0.8, seed=40
... )
... model.build()
... model.W
... 
array([[0.        , 0.12066975, 0.        , 0.        , 0.03891014],
       [0.12066975, 0.        , 0.        , 0.        , 0.98724568],
       [0.        , 0.10123466, 0.        , 0.17043792, 0.        ],
       [0.        , 0.01540023, 0.17043792, 0.        , 0.        ],
       [0.03891014, 0.98724568, 0.        , 0.        , 0.        ]])

Julia version:
1.

julia> rng = Random.seed!(42)
TaskLocalRNG()

julia> rand_hyper(rng, 5, 5; top_k=2)
5×5 Matrix{Float32}:
 0.0        0.0805158  0.0100425  0.0  0.0
 0.0805158  0.0        0.101827   0.0  0.0
 0.0        0.101827   0.0        0.0  0.990938
 0.0        0.0        0.0517202  0.0  0.0424772
 0.0        0.0734909  0.990938   0.0  0.0

julia> rng = Random.seed!(41)
TaskLocalRNG()

julia> rand_hyper(rng, 5, 5; top_k=2)
5×5 Matrix{Float32}:
 0.0  0.13659    0.0       0.0494492  0.0
 0.0  0.0        0.256201  0.363995   0.0
 0.0  0.256201   0.0       0.807059   0.0
 0.0  0.363995   0.807059  0.0        0.0
 0.0  0.0382711  0.0       0.0391867  0.0

julia> rng = Random.seed!(40)
TaskLocalRNG()

julia> rand_hyper(rng, 5, 5; top_k=2)
5×5 Matrix{Float32}:
 0.0       0.0  0.306278   0.0       0.0858886
 0.0       0.0  0.0821378  0.0       0.203433
 0.306278  0.0  0.0        0.0       0.756616
 0.0       0.0  0.20019    0.0       0.508322
 0.0       0.0  0.756616   0.508322  0.0