Piccolo matrix-free objective-constructor slice

src/control/objectives.jl

@@ -151,8 +151,54 @@ function CoherentKetInfidelityObjective(
         return abs(1 - coherent_ket_fidelity(ψ̃s, goals))
     end
 
+    # Declarable matrix-free per-knot HVP. The coherent loss is exactly
+    # `ℓ = |1 - F|` with `F = ‖A·z‖²`, so the per-knot Hessian factors as
+    # `Aᵀ · G · A` with `G = -2·sign(1 - F)·I_2` (the `:neg2_sign` rule).
+    A = _coherent_ket_lowrank_factor(goals, state_dims)
+    cap = ConstantLowRankHVP(A, :neg2_sign)
+
     # Pass vector of component names for multi-component terminal objective
-    return TerminalObjective(ℓ, ψ̃_names, traj; Q = Q)
+    return TerminalObjective(ℓ, ψ̃_names, traj; Q = Q, knot_hvp = cap)
+end
+
+# Build the constant rank-2 factor `A :: Matrix{Float64}` of size
+# `(2, sum(state_dims))` such that `[A·z]_1 = Re(S)`, `[A·z]_2 = Im(S)`,
+# where `S = (1/K) · Σᵢ ⟨gᵢ, ψᵢ⟩` is the coherent overlap and the
+# concatenated iso-state `z = [ψ̃₁; ψ̃₂; …; ψ̃_K]` follows the per-ket
+# `[real(ψᵢ); imag(ψᵢ)]` convention from `ket_to_iso`.
+# Then `‖A·z‖² = |S|² = F` exactly.
+function _coherent_ket_lowrank_factor(
+    goals::AbstractVector{<:AbstractVector{<:Complex}},
+    state_dims::AbstractVector{Int},
+)
+    K = length(goals)
+    @assert length(state_dims) == K
+    m = sum(state_dims)
+    A = zeros(Float64, 2, m)
+    inv_K = 1.0 / K
+    offset = 0
+    for i = 1:K
+        d_iso = state_dims[i]
+        @assert iseven(d_iso) "iso-state $i has odd dim $(d_iso); expected even"
+        d_c = d_iso ÷ 2
+        @assert length(goals[i]) == d_c (
+            "goal $i has $(length(goals[i])) entries; expected $(d_c) " *
+            "(half the iso-state dim)"
+        )
+        re_g = real(goals[i])
+        im_g = imag(goals[i])
+        # Block layout per ket: `ψ̃ᵢ = [Re(ψᵢ); Im(ψᵢ)]` (length d_iso).
+        #   ⟨gᵢ, ψᵢ⟩ = (Re(g) − i·Im(g))ᵀ (Re(ψ) + i·Im(ψ))
+        #             = Re(g)ᵀRe(ψ) + Im(g)ᵀIm(ψ)  +  i (Re(g)ᵀIm(ψ) − Im(g)ᵀRe(ψ))
+        # Row 1 (Re(S)) — scaled by 1/K
+        A[1, offset .+ (1:d_c)] .= inv_K .* re_g
+        A[1, offset .+ ((d_c+1):d_iso)] .= inv_K .* im_g
+        # Row 2 (Im(S)) — scaled by 1/K
+        A[2, offset .+ (1:d_c)] .= -inv_K .* im_g
+        A[2, offset .+ ((d_c+1):d_iso)] .= inv_K .* re_g
+        offset += d_iso
+    end
+    return A
 end
 
 # ---------------------------------------------------------
@@ -275,13 +321,52 @@ function unitary_fidelity_loss(Ũ⃗::AbstractVector{<:Real}, op::EmbeddedOpera
 end
 
 function UnitaryInfidelityObjective(
-    U_goal::AbstractPiccoloOperator,
+    U_goal::AbstractPiccoloOperator,  # full-op (AbstractMatrix) gets knot_hvp; EmbeddedOperator stays on fallback
     Ũ⃗_name::Symbol,
     traj::NamedTrajectory;
     Q = 100.0,
 )
+    # Declarable matrix-free per-knot HVP. The full-operator
+    # `ℓ = |1 − |tr(U†_goal · U)|²/n²|` factors exactly as `|1 − ‖A·Ũ⃗‖²|`
+    # with a constant rank-2 `A` (consumer-side rule `:neg2_sign`).
+    # The embedded-operator branch has a different (mixed-regularizer)
+    # form and stays on the dense fallback.
     ℓ = Ũ⃗ -> abs(1 - unitary_fidelity_loss(Ũ⃗, U_goal))
-    return TerminalObjective(ℓ, Ũ⃗_name, traj; Q = Q)
+    cap = if U_goal isa AbstractMatrix{<:Number}
+        ConstantLowRankHVP(_unitary_lowrank_factor(U_goal), :neg2_sign)
+    else
+        nothing
+    end
+    return TerminalObjective(ℓ, Ũ⃗_name, traj; Q = Q, knot_hvp = cap)
+end
+
+# Build the constant rank-2 factor `A :: Matrix{Float64}` of size
+# `(2, 2·n²)` such that `[A·Ũ⃗]_1 = Re(S)/n`, `[A·Ũ⃗]_2 = Im(S)/n`,
+# where `S = tr(U_goal† · U)` and the iso-vec `Ũ⃗` follows the
+# `operator_to_iso_vec` convention: per column i (0-indexed),
+# block of length 2n is `[Re(U[:, i+1]); Im(U[:, i+1])]`.
+# Then `‖A·Ũ⃗‖² = |S|²/n² = F` exactly.
+function _unitary_lowrank_factor(U_goal::AbstractMatrix{<:Number})
+    n = size(U_goal, 1)
+    @assert size(U_goal, 2) == n "U_goal must be square; got $(size(U_goal))"
+    A = zeros(Float64, 2, 2 * n^2)
+    inv_n = 1.0 / n
+    for i = 0:(n-1)
+        col = @view U_goal[:, i+1]
+        re_g = real(col)
+        im_g = imag(col)
+        # Per-column block of `Ũ⃗`: `i*2n + (1:n)` is Re, `i*2n + (n+1:2n)` is Im.
+        #   ⟨gᵢ, Uᵢ⟩ = (Re(g) − i·Im(g))ᵀ (Re(U) + i·Im(U))
+        #             = Re(g)ᵀRe(U) + Im(g)ᵀIm(U)  +  i (Re(g)ᵀIm(U) − Im(g)ᵀRe(U))
+        # Sum over columns gives S = Σᵢ ⟨gᵢ, Uᵢ⟩ = tr(U_goal† U).
+        # Row 1 (Re(S)/n)
+        A[1, i*2n .+ (1:n)] .= inv_n .* re_g
+        A[1, i*2n .+ ((n+1):2n)] .= inv_n .* im_g
+        # Row 2 (Im(S)/n)
+        A[2, i*2n .+ (1:n)] .= -inv_n .* im_g
+        A[2, i*2n .+ ((n+1):2n)] .= inv_n .* re_g
+    end
+    return A
 end
 
 function UnitaryFreePhaseInfidelityObjective(
@@ -424,11 +509,8 @@ using TestItems
     u = randn(1, N)
     Δt = fill(0.1, N)
 
-    traj = NamedTrajectory(
-        (ψ̃1 = ψ̃1, ψ̃2 = ψ̃2, u = u, Δt = Δt);
-        timestep = :Δt,
-        controls = :u,
-    )
+    traj =
+        NamedTrajectory((ψ̃1 = ψ̃1, ψ̃2 = ψ̃2, u = u, Δt = Δt); timestep = :Δt, controls = :u)
 
     # Goal states for X gate: |0⟩→|1⟩ and |1⟩→|0⟩
     ψ0 = ComplexF64[1.0, 0.0]
@@ -634,4 +716,204 @@ end
     gradient!(∇, obj, traj)
 end
 
+@testitem "knot_hvp matrix-free apply matches dense Hessian — CoherentKetInfidelityObjective" begin
+    using NamedTrajectories
+    using DirectTrajOpt
+    using DirectTrajOpt.Objectives: knot_hvp, ConstantLowRankHVP, get_full_hessian
+    using TrajectoryIndexingUtils
+    using LinearAlgebra
+    using Random
+
+    Random.seed!(42)
+    N = 10
+    ket_dim = 4  # iso dim for 2-level system
+    u = randn(1, N)
+    Δt = fill(0.1, N)
+    Q = 7.0  # non-default so the scaling threads through
+
+    # K=1 (single-ket) and K=2 (two-ket) — issue #257 AC requires ≥1 and ≥2 goals.
+    goal_bank = [ComplexF64[0.0, 1.0], ComplexF64[1.0, 0.0]]
+
+    for K in (1, 2)
+        goals_K = goal_bank[1:K]
+        names_K = K == 1 ? [:ψ̃1] : [:ψ̃1, :ψ̃2]
+
+        # Build a K-ket trajectory with a scaled copy of the goals at the
+        # terminal knot. `scale < 1` puts F below the kink; `scale > 1`
+        # puts F above.
+        function _make_traj(scale)
+            if K == 1
+                ψ̃1 = randn(ket_dim, N)
+                ψ̃1[:, N] .= scale .* ket_to_iso(goals_K[1])
+                return NamedTrajectory(
+                    (ψ̃1 = ψ̃1, u = u, Δt = Δt);
+                    timestep = :Δt,
+                    controls = :u,
+                )
+            else
+                ψ̃1 = randn(ket_dim, N)
+                ψ̃2 = randn(ket_dim, N)
+                ψ̃1[:, N] .= scale .* ket_to_iso(goals_K[1])
+                ψ̃2[:, N] .= scale .* ket_to_iso(goals_K[2])
+                return NamedTrajectory(
+                    (ψ̃1 = ψ̃1, ψ̃2 = ψ̃2, u = u, Δt = Δt);
+                    timestep = :Δt,
+                    controls = :u,
+                )
+            end
+        end
+
+        for scale in (0.4, 1.2)  # F < 1 and F > 1
+            traj = _make_traj(scale)
+            obj = CoherentKetInfidelityObjective(goals_K, names_K, traj; Q = Q)
+
+            cap = knot_hvp(obj, traj)
+            @test cap isa ConstantLowRankHVP
+            @test cap.core === :neg2_sign
+            @test size(cap.A) == (2, K * ket_dim)
+
+            # Index map for the terminal-knot concatenated iso-state.
+            knot_idx = vcat([slice(N, traj.components[nm], traj.dim) for nm in names_K]...)
+            z = vec(traj)
+            z_k = z[knot_idx]
+
+            # Matrix-free apply at the terminal knot:
+            #   H · v_k = Q · (-2·sign(1-F)) · Aᵀ · (A · v_k)
+            A = cap.A
+            F = sum(abs2, A * z_k)
+            G = -2.0 * sign(1.0 - F)
+            v = randn(length(z))
+            v_k = v[knot_idx]
+            Hv_mf = Q * G * (A' * (A * v_k))
+
+            # Dense apply at the same knot, via Symmetric(get_full_hessian, :U)·v.
+            H_full = Matrix(get_full_hessian(obj, traj))
+            Hv_dense = (Symmetric(H_full, :U)*v)[knot_idx]
+
+            @test Hv_mf ≈ Hv_dense atol = 1e-12
+        end
+    end
+end
+
+@testitem "knot_hvp matrix-free apply matches dense Hessian — UnitaryInfidelityObjective (full-op)" begin
+    using NamedTrajectories
+    using DirectTrajOpt
+    using DirectTrajOpt.Objectives: knot_hvp, ConstantLowRankHVP, get_full_hessian
+    using TrajectoryIndexingUtils
+    using LinearAlgebra
+    using Random
+
+    N = 10
+    u = randn(1, N)
+    Δt = fill(0.1, N)
+    Q = 5.0
+
+    # Issue #248 AC requires coverage across n = 2, 3, 4.
+    for n in (2, 3, 4)
+        Random.seed!(43 + n)             # per-n seed for reproducibility
+        iso_dim = 2 * n^2
+        # Random complex goal (no unitarity required for the parity test).
+        U_goal = randn(ComplexF64, n, n)
+
+        function _make_traj(scale)
+            Utilde = randn(iso_dim, N)
+            Utilde[:, N] .= scale .* operator_to_iso_vec(U_goal)
+            return NamedTrajectory(
+                (Utilde = Utilde, u = u, Δt = Δt);
+                timestep = :Δt,
+                controls = :u,
+            )
+        end
+
+        for scale in (0.4, 1.2)  # F < 1 and F > 1
+            traj = _make_traj(scale)
+            obj = UnitaryInfidelityObjective(U_goal, :Utilde, traj; Q = Q)
+
+            cap = knot_hvp(obj, traj)
+            @test cap isa ConstantLowRankHVP
+            @test cap.core === :neg2_sign
+            @test size(cap.A) == (2, iso_dim)
+
+            knot_idx = slice(N, traj.components[:Utilde], traj.dim)
+            z = vec(traj)
+            z_k = z[knot_idx]
+
+            A = cap.A
+            F = sum(abs2, A * z_k)
+            G = -2.0 * sign(1.0 - F)
+            v = randn(length(z))
+            v_k = v[knot_idx]
+            Hv_mf = Q * G * (A' * (A * v_k))
+
+            H_full = Matrix(get_full_hessian(obj, traj))
+            Hv_dense = (Symmetric(H_full, :U)*v)[knot_idx]
+
+            @test Hv_mf ≈ Hv_dense atol = 1e-12
+        end
+    end
+end
+
+@testitem "knot_hvp returns nothing for EmbeddedOperator UnitaryInfidelityObjective" begin
+    using NamedTrajectories
+    using DirectTrajOpt
+    using DirectTrajOpt.Objectives: knot_hvp
+    using LinearAlgebra
+
+    n_full = 4   # ambient
+    iso_dim = 2 * n_full^2
+    N = 5
+
+    Utilde = randn(iso_dim, N)
+    u = randn(1, N)
+    Δt = fill(0.1, N)
+    traj = NamedTrajectory((Utilde = Utilde, u = u, Δt = Δt); timestep = :Δt, controls = :u)
+
+    # 2x2 X gate embedded in a 4-dim ambient space (subspace = 1:2).
+    op = EmbeddedOperator(ComplexF64[0 1; 1 0], 1:2, n_full)
+    obj = UnitaryInfidelityObjective(op, :Utilde, traj; Q = 1.0)
+
+    # Embedded path stays on the dense fallback — no declared knot_hvp.
+    @test knot_hvp(obj, traj) === nothing
+end
+
+@testitem "knot_hvp field does not leak into KnotPointObjective value/gradient" begin
+    using NamedTrajectories
+    using DirectTrajOpt
+    using DirectTrajOpt.Objectives: ConstantLowRankHVP
+    using LinearAlgebra
+    using Random
+
+    # Issue #257 AC 3 / #248 AC 4: "Value and gradient of the objective are
+    # byte-for-byte unchanged — only the `knot_hvp` field is newly populated."
+    # Construct two `KnotPointObjective`s that differ ONLY in `knot_hvp`
+    # (one `nothing`, one a populated `ConstantLowRankHVP`) and assert
+    # `objective_value` and `gradient!` are literally equal on both.
+
+    Random.seed!(44)
+    N = 10
+    ket_dim = 4
+    ψ̃1 = randn(ket_dim, N)
+    u = randn(1, N)
+    Δt = fill(0.1, N)
+    traj = NamedTrajectory((ψ̃1 = ψ̃1, u = u, Δt = Δt); timestep = :Δt, controls = :u)
+
+    ℓ = x -> sum(x .^ 2) + 0.3 * sum(sin.(x))  # any non-trivial C² loss
+    Q = 7.0
+    fake_A = randn(2, ket_dim)                  # arbitrary; must not affect value/grad
+    cap = ConstantLowRankHVP(fake_A, :neg2_sign)
+
+    obj_none = KnotPointObjective(ℓ, :ψ̃1, traj; Qs = [Q], times = [N], knot_hvp = nothing)
+    obj_pop = KnotPointObjective(ℓ, :ψ̃1, traj; Qs = [Q], times = [N], knot_hvp = cap)
+
+    # Value: literally equal.
+    @test objective_value(obj_none, traj) === objective_value(obj_pop, traj)
+
+    # Gradient: literally equal (elementwise, no tolerance).
+    ∇_none = zeros(traj.dim * traj.N + traj.global_dim)
+    ∇_pop = zeros(traj.dim * traj.N + traj.global_dim)
+    gradient!(∇_none, obj_none, traj)
+    gradient!(∇_pop, obj_pop, traj)
+    @test ∇_none == ∇_pop
+end
+
 end

src/quantum/operators/direct_sums.jl

@@ -34,9 +34,7 @@ end
 Returns the direct sum of two iso_vec operators.
 """
 function direct_sum(Ã⃗::AbstractVector, B̃⃗::AbstractVector)
-    return operator_to_iso_vec(
-        direct_sum(iso_vec_to_operator(Ã⃗), iso_vec_to_operator(B̃⃗)),
-    )
+    return operator_to_iso_vec(direct_sum(iso_vec_to_operator(Ã⃗), iso_vec_to_operator(B̃⃗)))
 end
 
 """

src/quantum/primitives/isomorphisms.jl

@@ -76,8 +76,7 @@ function iso_vec_to_operator(Ũ⃗::AbstractVector{ℝ}) where {ℝ<:Real}
     N = Int(sqrt(Ũ⃗_dim))
     U = Matrix{complex(ℝ)}(undef, N, N)
     for i = 0:(N-1)
-        U[:, i+1] .=
-            @view(Ũ⃗[i*2N .+ (1:N)]) + one(ℝ) * im * @view(Ũ⃗[i*2N .+ ((N+1):2N)])
+        U[:, i+1] .= @view(Ũ⃗[i*2N .+ (1:N)]) + one(ℝ) * im * @view(Ũ⃗[i*2N .+ ((N+1):2N)])
     end
     return U
 end

src/quantum/systems/variational_quantum_systems.jl

@@ -22,11 +22,8 @@ The variational operators represent directions of uncertainty or perturbation in
 
 See also [`QuantumSystem`](@ref Piccolo.Quantum.QuantumSystems.QuantumSystem), [`OpenQuantumSystem`](@ref).
 """
-struct VariationalQuantumSystem{
-    F1<:Function,
-    F2<:Function,
-    F⃗3<:AbstractVector{<:Function},
-} <: AbstractQuantumSystem
+struct VariationalQuantumSystem{F1<:Function,F2<:Function,F⃗3<:AbstractVector{<:Function}} <:
+       AbstractQuantumSystem
     H::F1
     G::F2
     G_vars::F⃗3

src/quantum/templates/transmons/transmon_system.jl

@@ -440,8 +440,8 @@ function TransmonCavitySystem(;
 
     # Drift Hamiltonian
     H_drift =
-        Δ̃ * b̂' * b̂ - χ * â' * â * b̂' * b̂ - χ′ * b̂'^2 * b̂^2 * â' * â -
-        K_q * â'^2 * â^2 - K_c * b̂'^2 * b̂^2
+        Δ̃ * b̂' * b̂ - χ * â' * â * b̂' * b̂ - χ′ * b̂'^2 * b̂^2 * â' * â - K_q * â'^2 * â^2 -
+        K_c * b̂'^2 * b̂^2
 
     # Drive operators
     H_drives = [

test/test_utils.jl

@@ -88,12 +88,8 @@ function named_trajectory_type_1(; free_time = true)
             Δt = ([0.1], [0.30000000000000004]),
         )
     else
-        components = (
-            Ũ⃗ = data[1:8, :],
-            a = data[9:10, :],
-            da = data[11:12, :],
-            dda = data[13:14, :],
-        )
+        components =
+            (Ũ⃗ = data[1:8, :], a = data[9:10, :], da = data[11:12, :], dda = data[13:14, :])
         controls = (:dda,)
         timestep = 0.2
         bounds = (a = ([-1.0, -1.0], [1.0, 1.0]), dda = ([-1.0, -1.0], [1.0, 1.0]))