Sthosvd sketching

docs/src/tucker.md

@@ -23,6 +23,34 @@ core(tucker_res)
 factors(tucker_res)
 ```
 
+## Large tensors without explicit unfoldings
+
+For a fixed multilinear rank, ST-HOSVD can use an implicit randomized backend:
+
+```julia
+using Random
+
+tucker_res = tucker(
+    A,
+    (100, 40, 60);
+    method = :sthosvd,
+    svd_backend = :randomized,
+    processing_order = [2, 3, 1],
+    oversampling = 16,
+    power_iterations = 1,
+    block_columns = 65_536,
+    rng = MersenneTwister(0),
+)
+```
+Instead of constructing ``A_{(k)}``, the backend evaluates randomized projections
+and subspace iterations as tensor contractions. The Gaussian sketch is generated in 
+bounded blocks, and the input is not copied before its first mode projection. This
+substantially reduces memory when the tensor is large and the requested ranks are
+small relative to the mode dimensions.
+The exact backend remains the default and is generally preferable for small tensors,
+near-full ranks. Randomized results do not support `error_bound`, because discarded
+singular values are not computed.
+
 ## Tucker Docs
 
 ```@docs

src/api/tucker.jl

@@ -21,6 +21,33 @@ For `method = :hooi`:
     - `:sthosvd`: Uses ST-HOSVD to initialize the Tucker factors.
     - `TuckerResult`: Uses an existing Tucker decomposition as the initial point.
 
+For fixed-rank `method = :sthosvd`:
+
+* `svd_backend = :exact` preserves the deterministic dense implementation.
+* `svd_backend = :randomized` computes truncated mode subspaces through implicit
+  tensor contractions. It avoids materializing mode unfoldings and is intended for
+  large tensors whose target ranks are much smaller than their mode dimensions.
+* `oversampling = 16`, `power_iterations = 1`, `block_columns = 65_536`, and
+  `rng = Random.default_rng()` configure the randomized backend.
+
+```julia
+result = tucker(
+    A,
+    (100, 40, 60);
+    method = :sthosvd,
+    svd_backend = :randomized,
+    processing_order = [2, 3, 1],
+    oversampling = 16,
+    power_iterations = 1,
+    block_columns = 65_536,
+    rng = MersenneTwister(0),
+)
+```
+
+The randomized backend records empty per-mode singular-value vectors because it
+does not compute the complete discarded spectra. Consequently, [`error_bound`](@ref)
+is unavailable for those results; reconstruction-based error metrics remain valid.
+
 ## Example
 
 ```julia-repl

src/core/tensor_ops.jl

@@ -90,6 +90,161 @@ function mode_n_product(A::AbstractArray, U::AbstractMatrix, mode::Int)
     return permutedims(Bperm, invperm(perm))
 end
 
+# These helpers perform mode contractions from tensor index labels instead of
+# first constructing a mode unfolding. They are intentionally internal: the
+# public API remains `mode_n_product`, while scalable Tucker algorithms can use
+# these operations without materializing a tensor-sized `permutedims` copy.
+function _implicit_mode_product(
+    A::AbstractArray,
+    U::AbstractMatrix,
+    mode::Int;
+    block_columns::Int = 65_536,
+)
+    N = ndims(A)
+    1 <= mode <= N || throw(ArgumentError("mode must be in 1:$N; received $mode"))
+    block_columns >= 1 || throw(ArgumentError("block_columns must be positive"))
+    size(U, 2) == size(A, mode) || throw(
+        DimensionMismatch(
+            "size(U, 2)=$(size(U, 2)) must match size(A, mode)=$(size(A, mode))",
+        ),
+    )
+
+    tensor_labels = ntuple(identity, N)
+    new_label = N + 1
+    matrix_labels = (new_label, mode)
+    output_labels = ntuple(m -> m == mode ? new_label : m, N)
+    output_dims = ntuple(m -> m == mode ? size(U, 1) : size(A, m), N)
+    output = Array{promote_type(eltype(A), eltype(U))}(undef, output_dims)
+    block_lengths = _mode_block_lengths(size(A), mode, block_columns)
+    block_axes = ntuple(m -> _axis_chunks(axes(A, m), block_lengths[m]), N)
+
+    for block in Iterators.product(block_axes...)
+        block_indices = Tuple(block)
+        A_block = @view A[block_indices...]
+        partial = TensorOperations.tensorcontract(
+            output_labels,
+            U,
+            matrix_labels,
+            A_block,
+            tensor_labels,
+        )
+        output_indices = ntuple(m -> m == mode ? axes(output, m) : block_indices[m], N)
+        @views output[output_indices...] .= partial
+    end
+    return output
+end
+
+function _implicit_mode_cross(
+    A::AbstractArray,
+    B::AbstractArray,
+    mode::Int;
+    block_columns::Int = 65_536,
+)
+    N = ndims(A)
+    ndims(B) == N || throw(DimensionMismatch("A and B must have the same number of modes"))
+    1 <= mode <= N || throw(ArgumentError("mode must be in 1:$N; received $mode"))
+    block_columns >= 1 || throw(ArgumentError("block_columns must be positive"))
+    @inbounds for m = 1:N
+        m == mode && continue
+        size(A, m) == size(B, m) || throw(
+            DimensionMismatch(
+                "A and B must agree outside mode $mode; " *
+                "size(A, $m)=$(size(A, m)), size(B, $m)=$(size(B, m))",
+            ),
+        )
+    end
+
+    labels_A = ntuple(identity, N)
+    new_label = N + 1
+    labels_B = ntuple(m -> m == mode ? new_label : m, N)
+    output = zeros(promote_type(eltype(A), eltype(B)), size(A, mode), size(B, mode))
+    block_lengths = _mode_block_lengths(size(A), mode, block_columns)
+    block_axes = ntuple(m -> _axis_chunks(axes(A, m), block_lengths[m]), N)
+
+    for block in Iterators.product(block_axes...)
+        block_indices_A = Tuple(block)
+        block_indices_B = ntuple(m -> m == mode ? axes(B, m) : block_indices_A[m], N)
+        A_block = @view A[block_indices_A...]
+        B_block = @view B[block_indices_B...]
+        partial = TensorOperations.tensorcontract(
+            (mode, new_label),
+            A_block,
+            labels_A,
+            B_block,
+            labels_B,
+        )
+        output .+= partial
+    end
+    return output
+end
+
+@inline function _axis_chunks(axis::AbstractUnitRange, chunk_length::Int)
+    chunk_length >= 1 || throw(ArgumentError("chunk length must be positive"))
+    first_index = first(axis)
+    last_index = last(axis)
+    return [
+        start:min(start+chunk_length-1, last_index) for
+        start = first_index:chunk_length:last_index
+    ]
+end
+
+function _mode_block_lengths(dims::NTuple{N,Int}, mode::Int, block_columns::Int) where {N}
+    block_columns >= 1 || throw(ArgumentError("block_columns must be positive"))
+    lengths = ones(Int, N)
+    lengths[mode] = dims[mode]
+    remaining = block_columns
+    @inbounds for m = 1:N
+        m == mode && continue
+        lengths[m] = min(dims[m], max(remaining, 1))
+        remaining = max(div(remaining, lengths[m]), 1)
+    end
+    return Tuple(lengths)
+end
+
+"""
+    _implicit_mode_sketch(A, mode, sketch_rank, rng; block_columns)
+
+Compute `unfold_mode(A, mode) * Omega` for an implicit Gaussian matrix `Omega`
+without constructing either the mode unfolding or the complete random matrix.
+The non-mode columns are visited in tensor blocks containing at most roughly
+`block_columns` entries, and each partial random projection is evaluated as a
+tensor contraction.
+"""
+function _implicit_mode_sketch(
+    A::AbstractArray{T,N},
+    mode::Int,
+    sketch_rank::Int,
+    rng::AbstractRNG;
+    block_columns::Int,
+) where {T<:AbstractFloat,N}
+    1 <= mode <= N || throw(ArgumentError("mode must be in 1:$N; received $mode"))
+    sketch_rank >= 1 || throw(ArgumentError("sketch_rank must be positive"))
+
+    other_modes = [m for m = 1:N if m != mode]
+    block_lengths = _mode_block_lengths(size(A), mode, block_columns)
+    block_axes = ntuple(m -> _axis_chunks(axes(A, m), block_lengths[m]), N)
+    tensor_labels = ntuple(identity, N)
+    sketch_label = N + 1
+    omega_labels = Tuple(vcat(other_modes, sketch_label))
+    sketch = zeros(T, size(A, mode), sketch_rank)
+
+    for block in Iterators.product(block_axes...)
+        block_indices = Tuple(block)
+        A_block = @view A[block_indices...]
+        omega_dims = Tuple(vcat([size(A_block, m) for m in other_modes], sketch_rank))
+        omega = randn(rng, T, omega_dims)
+        partial = TensorOperations.tensorcontract(
+            (mode, sketch_label),
+            A_block,
+            tensor_labels,
+            omega,
+            omega_labels,
+        )
+        sketch .+= partial
+    end
+    return sketch
+end
+
 @inline function _rank1_entry_product(
     I::CartesianIndex{N},
     U::Vector{<:AbstractVector{T}},