Skip to contents

Project matrix onto a permutation-invariant space

Source: R/project_matrix.R
project_matrix.Rd

Project S onto the space of symmetric matrices invariant under the cyclic group generated by perm. Entries that must be equal under this symmetry are replaced by their average.

Usage

project_matrix(S, perm, precomputed_equal_indices = NULL)

Arguments

S

A square empirical covariance matrix to be projected. (See the S parameter in the gips() function). For optional multi-sample use, S can be a list of covariance matrices; see the Multi-sample section below. When it is not positive semi-definite, it shows a warning of a class not_positive_semi_definite_matrix.

perm

A permutation to be projected on. An object of a gips class, a gips_perm class, or anything that can be used as the x argument in the gips_perm() function.

precomputed_equal_indices

This parameter is no longer used and will be deleted in v1.4.0.

Value

The matrix S projected on the space of symmetric matrices invariant by a cyclic group generated by perm. See Covariance estimation for its statistical interpretation.

When S is a list, returns a list of projected matrices.

Details

The permutation can be specified in advance or selected using find_MAP(). Model selection by MAP concerns the cyclic group; the subsequent projection estimates the covariance matrix within the selected model.

The projection implements Definition 3 from references.

Covariance estimation

Let Z contain n observations in rows. The function returns the projection of the supplied S without changing its normalization:

  • If the mean is known to be zero and S = crossprod(Z) / n, the projection is the covariance MLE in the specified symmetry model, when that MLE exists.

  • If S = cov(Z), the projection retains the denominator n - 1 used by stats::cov(). For the full Gaussian likelihood maximized over both an unrestricted mean and a covariance matrix with the specified symmetry, the covariance MLE would be (n - 1) / n * project_matrix(S, perm), when it exists.

For S = cov(Z) with was_mean_estimated = TRUE, gips uses the projection without this rescaling, and logLik.gips() uses n - 1 in its likelihood calculation. This is the package's convention based on centered data, distinct from the full Gaussian likelihood described above. See the Estimated mean subsection in vignette("Theory", package = "gips").

MLE availability

project_matrix() does not receive the sample size and does not check whether an MLE is available. It can return a singular projected matrix. For a corresponding gips object g, summary.gips() reports the sample-size threshold in summary(g)$n0. The field summary(g)$MLE_available combines the sample-size check with a numerical positive-definiteness check of the projection, using the package's convention described above.

Multi-sample

When S is a list of matrices (as stored in a multi-sample gips object), project_matrix() applies the same permutation projection to each element and returns a list of projected matrices of the same length.

See also

  • Wikipedia - Estimation of covariance matrices

  • Project Matrix - Equation (6) section of vignette("Theory", package = "gips") or its pkgdown page - A place to learn more about the math behind the gips package and see more examples of project_matrix().

  • find_MAP() - Searches for a cyclic group with maximum posterior weight. Heuristic optimizers do not guarantee finding the maximum. The returned permutation can be used to project S onto the selected symmetry model (see examples).

  • gips_perm() - Constructor for the perm parameter.

  • plot.gips() - For plot(g, type = "MLE"), the project_matrix() is called (see examples).

  • summary.gips() - Reports n0 and MLE_available for a gips object; see the MLE availability section.

Examples

p <- 6
my_perm <- "(14)(23)" # permutation (1,4)(2,3)(5)(6)
number_of_observations <- 10
X <- matrix(rnorm(p * number_of_observations), number_of_observations, p)
S <- cov(X)
projected_S <- project_matrix(S, perm = my_perm)
projected_S
#>             [,1]        [,2]        [,3]        [,4]        [,5]       [,6]
#> [1,]  1.09598470 -0.06413469 -0.21185653 -0.43353246  0.05920705  0.2202124
#> [2,] -0.06413469  1.00959419  0.09589337 -0.21185653  0.48366779 -0.3530059
#> [3,] -0.21185653  0.09589337  1.00959419 -0.06413469  0.48366779 -0.3530059
#> [4,] -0.43353246 -0.21185653 -0.06413469  1.09598470  0.05920705  0.2202124
#> [5,]  0.05920705  0.48366779  0.48366779  0.05920705  1.15745710 -0.2385771
#> [6,]  0.22021240 -0.35300588 -0.35300588  0.22021240 -0.23857709  0.7438371
# The value in [1,1] is the same as in [4,4]; also, [2,2] and [3,3];
  # also [1,2] and [3,4]; also, [1,5] and [4,5]; and so on

# Plot the projected matrix:
g <- gips(S, number_of_observations, perm = my_perm)
plot(g, type = "MLE")


# Search for the MAP symmetry group:
g_MAP <- find_MAP(g, max_iter = 10, show_progress_bar = FALSE, optimizer = "Metropolis_Hastings")
# Then, project S onto the selected model:
S_MAP <- project_matrix(attr(g, "S"), perm = g_MAP)
S_MAP
#>              [,1]        [,2]         [,3]         [,4]        [,5]
#> [1,]  0.686203184  0.06033489 -0.250237627 -0.009842726  0.06033489
#> [2,]  0.060334885  1.28496383  0.060334885 -0.317893688  0.74598323
#> [3,] -0.250237627  0.06033489  0.686203184 -0.009842726  0.06033489
#> [4,] -0.009842726 -0.31789369 -0.009842726  1.483914734 -0.31789369
#> [5,]  0.060334885  0.74598323  0.060334885 -0.317893688  1.28496383
#> [6,] -0.250237627  0.06033489 -0.250237627 -0.009842726  0.06033489
#>              [,6]
#> [1,] -0.250237627
#> [2,]  0.060334885
#> [3,] -0.250237627
#> [4,] -0.009842726
#> [5,]  0.060334885
#> [6,]  0.686203184
plot(g_MAP, type = "heatmap")