Function Repository Resource:

CenteringMatrix

Source Notebook

Generate an n-sized centering matrix

Contributed by: Phileas Dazeley-Gaist

ResourceFunction["CenteringMatrix"][n]

generates an n-sized centering matrix.

Details and Options

The centering matrix of size n is defined as the n×n matrix:

, where I_n is the identity matrix of size n and J_n is an n×n matrix of ones.

Centering a vector means to subtract the mean of the vector from every one of its element.

Column-wise centering of a matrix M means subtracting the mean of each column from its elements. This is achieved by left-multiplication of the centering matrix C by M, that is, forming C·M.

Column-wise centering of a matrix M means subtracting the mean of each row from its elements. This is achieved by right-multiplication of the centering matrix C by M, M·C.

By default, the function returns a SparseArray of positive machine precision real values. Specify the working precision with the WorkingPrecision option (defaults to MachinePrecision).

When centering large vectors or matrices, it is good practice to use numerical approximation, either by setting the WorkingPrecision, or by applying N to the matrix being processed.

Examples

Basic Examples (5)

Generate the 3×3 centering matrix:

In[1]:=

Out[1]=

Center a vector:

In[2]:=

Out[3]=

Confirm the vector is centered:

In[4]:=

Out[4]=

Center a matrix column-wise:

In[5]:=

m = Array[RandomInteger[{1, 10}] &, {3, 5}];
mCentered = ResourceFunction["CenteringMatrix"][Length[m]] . m;
MatrixForm[mCentered]

Out[4]=

Confirm the matrix is centered:

In[6]:=

Out[6]=

Center a matrix row-wise:

In[7]:=

m = Array[RandomInteger[{1, 10}] &, {3, 5}]; mCentered = m . ResourceFunction["CenteringMatrix"][Last[Dimensions[m]]];
MatrixForm[mCentered]

Out[8]=

Confirm the matrix is centered:

In[9]:=

Out[9]=

Double center a matrix:

In[10]:=

m = Array[RandomInteger[{1, 10}] &, {3, 5}]; mCentered = ResourceFunction["CenteringMatrix"][Length[m]] . m . ResourceFunction["CenteringMatrix"][Last[Dimensions[m]]];
MatrixForm[mCentered]

Out[11]=

Confirm the matrix is double centered:

In[12]:=

$N@Round[#, 10^-10] &@{Mean[mCentered], Mean[Transpose[mCentered]]}$

Out[12]=

Options (6)

WorkingPrecision (4)

By default WorkingPrecision is set to MachinePrecision:

In[13]:=

Out[13]=

Compute the 3×3 centering matrix with 2-digit precision:

In[14]:=

Out[14]=

Compute the 3×3 centering matrix with machine precision:

In[15]:=

Out[15]=

Compute the 3×3 centering matrix with exact values:

In[16]:=

$ResourceFunction["CenteringMatrix"][3, WorkingPrecision -> \[Infinity]] // MatrixForm$

Out[16]=

TargetStructure (2)

By default, CenteringMatrix will return a SparseArray:

In[17]:=

Out[17]=

Explicitly request a dense matrix from CenteringMatrix:

In[18]:=

Out[18]=

Applications (3)

Define a fictional tabular dataset of species observations by site:

In[19]:=

dataset = Transpose@
Dataset@AssociationThread[
Table["Species " <> ToUpperCase@Alphabet[][[i]], {i, 4}] -> Ceiling@RandomVariate[NormalDistribution[2, 1], {4, 6}]]

Out[19]=

Double center the dataset:

In[20]:=

datasetCentered = Module[{keys = Normal[dataset[First, Keys]], values = Normal[Values[dataset]], dims},
dims = Dimensions[values];
values = ResourceFunction["CenteringMatrix"][First[dims]] . values . ResourceFunction["CenteringMatrix"][Last[dims]];
Transpose@Dataset@AssociationThread[keys -> Transpose[values]]
]

Out[20]=

Confirm the dataset is double centered:

In[21]:=

$With[{vals = Values[datasetCentered]}, Round[#, 10^-10] &@{Mean[vals], Mean[Transpose[vals]]}]$

Out[21]=

Neat Examples (4)

Define a function to perform PCoA (Classical Multidimensional Scaling) that uses CenteringMatrix:

In[22]:=

$ClearAll[PCoA] PCoA[distanceMatrix_List, dim : _Integer ? Positive : 2] := Module[{centered, c, eigVal, eigVect}, (*Centering matrix*) c = ResourceFunction["CenteringMatrix"][Length[distanceMatrix], WorkingPrecision -> MachinePrecision]; (*Double centered distance matrix*) centered = -(1/2) c . distanceMatrix . c; (*Eigenvalues and eigenvectors*){eigVal, eigVect} = Eigensystem[centered, dim]; (*Coordinates*) Transpose[eigVect] . Sqrt[DiagonalMatrix[eigVal]] // Re ]$