Solve Weighted Procrustes Rotation Problem — solve_procrustes_rotation

Finds the optimal rotation matrix `R` that aligns a source matrix `A_source` to a target matrix `T_target`, minimizing `|| Omega (A_source R - T_target) ||_F^2`, where `Omega` is a diagonal weighting matrix.

Usage

solve_procrustes_rotation_weighted(
  A_source,
  T_target,
  m_parcel_rows,
  m_task_rows,
  omega_mode = "fixed",
  fixed_omega_weights = NULL,
  reliability_scores = NULL,
  scale_omega_trace = TRUE
)

Arguments

A_source: The numeric source matrix (e.g., augmented anchors for one subject), dimensions `N x k` (N = total anchors, k = dimensions).
T_target: The numeric target matrix (e.g., augmented template), dimensions `N x k`.
m_parcel_rows: Integer, number of rows in `A_source` (and `T_target`) corresponding to parcel anchors. Assumed to be the first `m_parcel_rows`.
m_task_rows: Integer, number of rows in `A_source` (and `T_target`) corresponding to task condition anchors. Assumed to be after parcel rows.
omega_mode: Character string, mode for determining weights: `"fixed"` (default) or `"adaptive"`.
fixed_omega_weights: List, used if `omega_mode == "fixed"`. E.g., `list(parcel = 1.0, condition = 0.5)`. If `NULL`, defaults to this.
reliability_scores: Numeric vector of length `m_task_rows`, containing reliability scores (e.g., R^2_p) for each task condition. Required and used only if `omega_mode == "adaptive"`.
scale_omega_trace: Logical, whether to rescale the diagonal `Omega` matrix so that `sum(diag(Omega)) == N`. Default `TRUE`.

Value

A `k x k` rotation matrix `R`.

Details

Internally computes the weighted cross-product using `crossprod(A_source, T_target * omega_diag_vector)` to avoid creating separate weighted matrices for `A_source` and `T_target`.

When the spectral dimension `k` is 1, the SVD-based determinant correction is unnecessary. In this case the rotation is simply the sign of the scalar cross-product `M`, i.e. `R <- matrix(sign(M), 1, 1)`.

Examples

N_parcels <- 5
N_tasks <- 3
N_total <- N_parcels + N_tasks
k_dims <- 4
A <- matrix(rnorm(N_total * k_dims), N_total, k_dims)
T_true <- matrix(rnorm(N_total * k_dims), N_total, k_dims)
R_true <- svd(matrix(rnorm(k_dims*k_dims),k_dims,k_dims))$u %*% 
          diag(sample(c(1,1,1,-1),k_dims,replace=TRUE)) %*% 
          svd(matrix(rnorm(k_dims*k_dims),k_dims,k_dims))$v
A_rotated_perfect <- T_true %*% t(R_true) # A_rotated_perfect R_true should be T_true

# Fixed weights (default)
R_fixed <- solve_procrustes_rotation_weighted(A_rotated_perfect, T_true, 
                                              N_parcels, N_tasks)
# print(all.equal(T_true %*% t(R_fixed), A_rotated_perfect)) # Check alignment
print(sum((A_rotated_perfect %*% R_fixed - T_true)^2)) # Should be small
#> [1] 9.915003

# Adaptive weights (example with dummy reliabilities)
set.seed(123)
rel_scores <- runif(N_tasks, -0.2, 0.8)
R_adaptive <- solve_procrustes_rotation_weighted(A_rotated_perfect, T_true, 
                                                 N_parcels, N_tasks, 
                                                 omega_mode = "adaptive", 
                                                 reliability_scores = rel_scores)
print(sum((A_rotated_perfect %*% R_adaptive - T_true)^2)) # Should be small
#> [1] 17.43792

# No task rows
R_no_task <- solve_procrustes_rotation_weighted(A[1:N_parcels,], T_true[1:N_parcels,],
                                                N_parcels, 0)