VLSVT

Size and shape of multivariate observations

Let \({\mathbf{Y}} = [y_{ij}]_{n \times q}\) be a matrix containing the values of \(q\) non-negative variables observed in \(n\) individuals.

# load("toy.RData")
# Y = t(data[c(1,2), -c(1,2)]) # Get only 1 gene
# covariate = metadata[, 3]

# hardcoded
Y <- cbind(
  c(4.66, 0.26, 3.46, 13.36, 0.87, 0.35, 14.29, 11.87, 16.68, 37.34, 7.59, 1.74, 10.17, 3.58, 20.39, 12.43, 10.64, 15.79),
  c(0.92, 0.88, 3.74, 5.37, 0.38, 0.39, 9.05, 5.48, 10.04, 3.25, 2.42, 0.00, 6.76, 3.61, 4.75, 3.29, 3.88, 4.58)
  )
covariate <- factor(c(rep("A", 6), rep("B", 6), rep("C", 6)))
  
n <- nrow(Y)
q <- ncol(Y)
rownames(Y) <- NULL
colnames(Y) <- paste0("y", 1:2)

y1	y2
4.66	0.92
0.26	0.88
3.46	3.74
13.36	5.37
0.87	0.38
0.35	0.39
14.29	9.05
11.87	5.48
16.68	10.04
37.34	3.25
7.59	2.42
1.74	0.00
10.17	6.76
3.58	3.61
20.39	4.75
12.43	3.29
10.64	3.88
15.79	4.58

We define size and shape as the overall and relative levels of the multivariate observations. For a given individual \(i\):

• Size: \[\lambda_i = \sum_{j=1}^q y_{ij}\]
• Shape: \[\mathbf{p}_i = \frac{\mathbf{y}_i}{\lambda_i}\] where \(\mathbf{y}_i\) is the vector of values of the \(q\) variables for individual \(i\).

For example, consider neuroimaging endophenotypes such as the MRI-derived volumes of brain regions. In this example, the size would be the volume of the whole brain, and the shape, the relative volume of the brain regions. In another example, consider the expression of the transcript isoforms of a given gene. The size would be the overall expression of the gene, and the shape the relative abundances of its transcripts.

A model for constant shape

We define the total variability of \(\mathbf{Y}\) (\(V_t\)) as the sum of the variances of the \(q\) variables across the \(n\) individuals (equivalently, the trace of the covariance matrix of the \(q\) variables):

\[ V_t = tr(\mathbf{\Sigma}) = \sum_{j=1}^q\sigma^2_j \tag{1}\]

where \({\mathbf{\Sigma}} = [\sigma_{jk}]_{q\times q}\) is the sample covariance matrix of the \(q\) variables. We may also work in terms of sums of squares: \(SS_t = V_t(n-1)\).

sigma <- cov(Y) * (n - 1) # crossprod(scale(Y, center = T, scale = F))
ss.t <- sum(diag(sigma)) # Working in terms of SS instead of variances

We define \(V_{ls}\) as the variability over the line of \(\mathbb{R}^q\) with non-negative coefficients that minimizes the distance between the original and projected \(n\) points. This is the line which better fits the values of the \(q\) variables observed in the \(n\) individuals under the following multiplicative model [@Gonzalez2012].

Our multiplicative model specifies that the value of a variable for a given individual (\(y_{ij}\)) is the product of the size parameter for this individual (\(\lambda_i\)) and the \(j\)th component of a shape parameter vector (\(\mathbf{p}\)) that is assumed constant for the \(n\) individuals:

\[ \mathbf{y}_i = \begin{pmatrix} y_{i1} \\ y_{i2} \\ \vdots \\ y_{iq} \end{pmatrix} = \lambda_i \begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_q \end{pmatrix} = \lambda_i \mathbf{p} \tag{2}\]

The ratio \(V_{ls}/V_t\) represents the fraction of total variability explained by a constant shape model.

Calculation of \(V_{ls}\)

We use the least squares criterion to estimate \(\mathbf{p}\) and \(\lambda_i\), restricted to \(p_j \ge 0\) and \(\lambda_i \ge 0\). Note that \(\mathbf{Y}^\textrm{T}\) here is \(\mathbf{Y}\) in @Gonzalez2012.

\[ \min_{\boldsymbol{\lambda} \geq 0,\,\\ \mathbf{p} \geq 0} \left\| \mathbf{Y} - \boldsymbol{\lambda} \mathbf{p}^\textrm{T} \right\|_F^2 \] Where \(\boldsymbol{\lambda} = [\lambda_j]_{n \times 1}\) is the vector of size parameters, and \(\| \cdot \|_F\) is the Frobenius norm.

Let \(d_1\), \(\mathbf{u}_1\) and \(\mathbf{v}_1\) denote the greater singular value of \(\mathbf{Y}\) and its corresponding left and right singular vectors. Given that \(\mathbf{Y}\) is a non-negative rectangular matrix, the Perron-Frobenius theorem (@Gantmacher1959) ensures that \(d_1\) and all the components \(u_{1i}\) of \(\mathbf{u}_1\) and \(v_{1j}\) of \(\mathbf{v}_1\) are non-negative. The Eckart-Young theorem (@EckartYoung1936) ensures that:

\[ \min_{\boldsymbol{\lambda} \geq 0,\,\\ \mathbf{p} \geq 0} \left\| \mathbf{Y} - \boldsymbol{\lambda} \mathbf{p}^\textrm{T} \right\|_F = \left\| \mathbf{Y} - d_1\mathbf{u}_1\mathbf{v}_1^\textrm{T} \right\|_F \]

The explicit expressions of \(\boldsymbol{\lambda}\) and \(\mathbf{p}\) are:

\[ \boldsymbol{\lambda} = \left( \sum_{j=1}^{q} v_{1j} \right) d_1 \mathbf{u}_1 \quad \text{and} \quad \mathbf{p} = \left( \sum_{j=1}^{q} v_{1j} \right)^{-1} \mathbf{v}_1 \]

d.svd <- svd(Y) # Check if u1 (v1) components have same sign
d1 <- d.svd$d[1]
u1 <- d.svd$u[, 1, drop = F]
v1 <- d.svd$v[, 1, drop = F]
lambda <- sum(v1) * d1 * u1 # sizes
p <- v1/sum(v1) # shape

The computation of \(V_{ls}\) requires only the values of \(v_1\) (the first right singular vector):

\[V_{ls} = \mathbf{v}_1^\textrm{T}\mathbf{\Sigma}\mathbf{v}_1 \tag{3}\]

Z <- lambda %*% t(p) # u1 %*% d1 %*% t(v1)
colnames(Z) <- colnames(Y)
rownames(Z) <- rownames(Y)
sigma.z <- cov(Z)
sigma.z <- sigma.z * (n - 1)
ss.t.z <- sum(diag(sigma.z)) # t(v1)%*%sigma%*%v1

The model of constant shape projects the original observations in the subspace of dimension 1 defined by the first right singular vector \(\mathbf{v}_1\). The projected data is \(\mathbf{Z} = \mathbf{Yv}_1\mathbf{v}_1^\textrm{T} = \mathbf{YP}\).

P <- tcrossprod(v1) # crossprod(v1) = 1; sqrt(sum(v1^2)) = 1 (orthonormal)
# all.equal(Y %*% P, Z) 

\(\mathbf{P}\) is idempotent (\(\mathbf{P}^2 = \mathbf{P}\)) and symmetric (\(\mathbf{P}^\textrm{T} = \mathbf{P}\)), and hence it defines an orthogonal projection.

Proofs

Given the singular value decomposition of \(\mathbf{Y}\): \[ \mathbf{Y} = \mathbf{UDV}^\textrm{T} = \sum_{k=1}^rd_k\mathbf{u}_k\mathbf{v}_k^\textrm{T} \] where \(r = rank(\mathbf{Y})\).

\[ \mathbf{Yv}_1 = \sum_{k=1}^rd_k\mathbf{u}_k\mathbf{v}_k^\textrm{T}\mathbf{v}_1 = d_1\mathbf{u}_1 \] given that \(\mathbf{v}_k^\textrm{T}\mathbf{v}_l = 0 \; \forall \; l \ne k\) and \(\mathbf{v}_k^\textrm{T}\mathbf{v}_k = 1\), as the \(\mathbf{v}_k\) are orthonormal.

Therefore \[\mathbf{Yv}_1\mathbf{v}_1^\textrm{T} = d_1\mathbf{u}_1\mathbf{v}_1^\textrm{T} = \mathbf{Z} \tag{4}\] If \(\mathbf{Z}\) is column centered (if it is not the case, the proof holds as well with some modifications),

\[tr(\boldsymbol{\Sigma}_{\mathbf{Z}}) = \frac{1}{n-1}\mathbf{Z}^\textrm{T}\mathbf{Z}\]

Substituting,

\[ \begin{align} tr(\boldsymbol{\Sigma}_{\mathbf{Z}}) = tr\left(\frac{1}{n-1}\left(\mathbf{YP}\right)^\textrm{T}\left(\mathbf{YP}\right)\right) = tr\left(\frac{1}{n-1}\mathbf{P}^\textrm{T}\mathbf{Y}^\textrm{T}\mathbf{Y}\mathbf{P}\right) = tr\left(\mathbf{P}^\textrm{T}\boldsymbol{\Sigma}\mathbf{P}\right) = \\ tr\left(\boldsymbol{\Sigma}\mathbf{P}\mathbf{P}^\textrm{T}\right) = tr\left(\boldsymbol{\Sigma}\mathbf{v}_1\mathbf{v}_1^\textrm{T}\mathbf{v}_1\mathbf{v}_1^\textrm{T}\right) = tr\left(\boldsymbol{\Sigma}\mathbf{v}_1\mathbf{v}_1^\textrm{T}\right) = tr\left(\mathbf{v}_1^\textrm{T}\boldsymbol{\Sigma}\mathbf{v}_1\right) = \mathbf{v}_1^\textrm{T}\boldsymbol{\Sigma}\mathbf{v}_1 \end{align} \]

\(V_{ls}/V_t =\) 0.9212682.

Covariates

The classical ANOVA decomposition of the total variability in its between groups and within groups components \[ SS_t = SS_b + SS_w \] can be generalized to the multivariate scenario based on the direct sum of the univariate sums of squares across all the variables [@Anderson2001].

From Equation 1, we have

\[ SS_t = (n-1)tr(\boldsymbol{\Sigma}) \]

In addition, it can be shown that

\[ SS_w = \sum_{k=1}^K(n_k-1) tr(\boldsymbol{\Sigma}_k) \]

Where \(n_k\) and \(\Sigma_k\) are, respectively, the sample size and sample covariance matrix for group \(k\) out of \(K\).

Hence,

\[ SS_b = (n-1)tr(\boldsymbol{\Sigma}) - \sum_{k=1}^K(n_k-1) tr(\boldsymbol{\Sigma}_k) \]

Considering again the constant shape model and the variability of the projected data (Equation 3), if \(SS_{b(Z)}\) is the corresponding between variability of the projected points (see also @Gonzalez2012SM), then

\[ \begin{align} SS_{b(Z)} = (n-1)tr(\boldsymbol{\Sigma}_{\mathbf{Z}}) - \sum_{k=1}^K(n_k-1) tr(\boldsymbol{\Sigma}_{\mathbf{Z}_k}) = \\ (n-1)\mathbf{v}_1^\textrm{T}\boldsymbol{\Sigma}\mathbf{v}_1 - \sum_{k=1}^K(n_k-1)\mathbf{v}_1^\textrm{T}(\boldsymbol{\Sigma}_k)\mathbf{v}_1 \end{align} \]

dz <- d <- 0
groups <- levels(covariate)
for(gp in groups){
  sel <- covariate == gp
  sigma.g <- cov(Y[sel, ]) * (sum(sel) - 1)
  d <- d + sum(diag(sigma.g))
  dz <- dz + t(v1)%*%sigma.g%*%v1
}

ss.b <- ss.t - d
ss.b.z <- ss.t.z - dz

# ss.b/ss.t
# ss.b.z/ss.b

The ratios of interest are \(\frac{SS_b}{SSt} =\) 0.275994 (the fraction of total variability explained by the covariate) and \(\frac{SS_b(\mathbf{Z})}{SS_b} =\) 0.999655 (the fraction of the variability explained by the covariate that is explained by a constant shape model).

Equivalently, from a linear modeling perspective, \[ SS_b = tr\left((\mathbf{H}\tilde{\mathbf{Y}})^\textrm{T}(\mathbf{H}\tilde{\mathbf{Y}})\right) = tr\left(\tilde{\mathbf{Y}}^\textrm{T}\mathbf{H}\tilde{\mathbf{Y}} \right) \] where \(\tilde{\mathbf{Y}}\) is column-centered \(\mathbf{Y}\) and \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\textrm{T}\mathbf{X})^{-1}\mathbf{X}^\textrm{T}\) is the well-known hat matrix. \(\mathbf{X}\) is the design matrix. Note that \(\mathbf{\hat{Y}}=\mathbf{H}\tilde{\mathbf{Y}}\) are the fitted values.

lm.fit <- lm(Y ~ covariate)
Y.lm <- lm.fit$fitted.values 
ss.b <- sum(diag(cov(Y.lm))) * (n - 1)  
Y.lm.vls <- Y.lm %*% P
ss.b.z <- sum(diag(cov(Y.lm.vls))) * (n - 1) 

# ss.b/ss.t
# ss.b.z/ss.b

# # Benchmark running time
# Zc <- scale(Z, center = T, scale = F)
# M <- model.matrix(~0+covariate)
# H <- M%*%solve(crossprod(M))%*%t(M)
# sum(diag(t(Zc)%*%H%*%Zc)) # ss.b.z
# sum(diag(cov(H%*%Z)))*(n-1) ## also; HZ = fitted values

# Yc <- scale(Y, center = T, scale = F)
# sum(diag(t(Yc)%*%H%*%Yc)) ## ss.b
# sum(diag(cov(H%*%Y)))*(n-1) ## also; HY = fitted values

This enables to incorporate continuous covariates. Note that the order of the projections is relevant: the meaning of the ratios changes. E.g. is not the same to compute the fraction of the variance explained by a covariate that can be explained by a constant shape model, than the fraction of variance explained by a constant shape model that can be explained by a covariate.

lm.fit <- lm(Z ~ covariate)
Y.vls.lm <- lm.fit$fitted.values 
ss.b.z <- sum(diag(cov(Y.vls.lm))) * (n - 1)

# ss.t.z/ss.t # vls
# ss.b.z/ss.t.z

A model for constant size

In contrast to the \(V_{ls}\) model (Equation 2), let’s consider now a model where the value of a variable for a given individual (\(y_{ij}\)) is the product of a size parameter that is assumed constant for all the individuals (\(\lambda\)), and the \(j\)th component of a shape parameter vector corresponding to this individual (\(\mathbf{p}_i\)).

\[ \mathbf{y}_i = \begin{pmatrix} y_{i1} \\ y_{i2} \\ \vdots \\ y_{iq} \end{pmatrix} = \lambda \begin{pmatrix} p_{i1} \\ p_{i2} \\ \vdots \\ p_{iq} \end{pmatrix} = \lambda \mathbf{p}_i \] If the parameters are renamed \(\lambda_{ij} = \lambda p_{ij}\), the above model becomes a linear model with constraints.

\[ y_{ij} = \lambda_{ij} \quad \text{with} \quad \sum_{j=1}^{q} \lambda_{1j} = %\sum_{j=1}^{q} \lambda_{2j} = \cdots = \sum_{j=1}^{q} \lambda_{nj} \quad \text{and} \quad \lambda_{ij} \ge 0 \tag{5}\]

This model can be interpreted geometrically as follows. The individuals are projected on a \((q-1)-\)simplex. Hence, the projections are non-negative and have a constant sum of components (i.e. constant size). The position of the simplex can be estimated by the least squares method. We define \(V_{ls(2)}\) as the variability over the \((q-1)-\)simplex with non-negative values that minimizes the distance between the original and projected \(n\) points.

A brief note about the simplex

The standard \((q-1)-\)simplex is a subset of \(\mathbb{R}^q\) that fulfils: \[ \Delta^{q-1} = \left\{ \mathbf{x} \in \mathbb{R}^q \;\Big|\; x_j \ge 0 \;\forall j = 1, \dots, q ,\;\sum_{j=1}^{q} x_j = 1 \right\} \] Here, we refer to a scaled \((q-1)-\)simplex of size \(\lambda > 0\), given by:

\[ \Delta^{q-1}_s = \left\{ \mathbf{x} \in \mathbb{R}^q \;\Big|\; x_j \ge 0 \;\forall j = 1, \dots, q ,\;\sum_{j=1}^{q} x_j = \lambda \right\} \tag{6}\]

Calculation of \(V_{ls(2)}\) (I)

Let’s not consider for now the \(np\) inequality constraints in Equation 5. We will minimize the function

\[ \mathcal{L}(\lambda_{ij}, \mu_i) = \sum_{i=1}^n \sum_{j=1}^q \bigl(y_{ij} - \lambda_{ij}\bigr)^2 + \sum_{i=1}^n \mu_i \left( \sum_{j=1}^q \lambda_{ij} - \lambda \right). \] with the constraint

\[ \sum_{j=1}^{q} \lambda_{ij} = \lambda \;\; \forall \; i = 1, \ldots, n \tag{7}\]

using a Lagrange multiplier. We obtain

\[ \lambda_{ij} = y_{ij} + \frac{1}{nq}\sum_{r=1}^n \sum_{s=1}^q y_{rs} - \frac{1}{q}\sum_{s=1}^q y_{is} \]

Proofs

Compute the partial derivative with respect to \(\lambda_{ij}\) and set it to 0: \[ \begin{align} \frac{\partial \mathcal{L}}{\partial \lambda_{ij}} = -2(y_{ij} - \lambda_{ij}) + \mu_i = 0\\ \Rightarrow \lambda_{ij} = y_{ij} + \frac{\mu_i}{2} \end{align} \]

Substitute \(\lambda_{ij}\) on the constraint equation (Equation 7):

\[ \sum_{j=1}^q \lambda_{ij} = \sum_{j=1}^q ( y_{ij} + \frac{\mu_i}{2} ) = \sum_{j=1}^q y_{ij} + q \cdot \frac{\mu_i}{2} = \lambda \]

Solve for \(\mu_i\):

\[ \mu_i = \frac{2}{q} \left( \lambda - \sum_{j=1}^q y_{ij} \right) \]

Substitute \(\mu_i\) in the first equation:

\[ \lambda_{ij} = y_{ij} + \frac{1}{q} \left ( \lambda - \sum_{j=1}^{q} y_{ij} \right) \]

We need to assume that \(\boldsymbol{\lambda}\) is the average size value across individuals:

\[ \lambda = \frac{1}{n} \sum_{i=1}^{n} \sum_{j=1}^{q} y_{ij} \] Hence

\[ \lambda_{ij} = y_{ij} + \frac{1}{nq} \sum_{i=1}^{n} \sum_{j=1}^{q} y_{ij} - \frac{1}{q} \sum_{j=1}^{q} y_{ij} \]

Let \(\bar{y}_i = \frac{1}{q}\sum_{j=1}^qy_{ij}\) be the average of a given row. Therefore

\[ \lambda_{ij} = y_{ij} + \left( \frac{\lambda}{q} - \bar{y}_i \right) \]

Interpretation: Each \(\lambda_{ij}\) is \(y_{ij}\) adjusted by a constant shift, given by the difference between the overall mean and the mean of the \(i\)th row. The \(\lambda_{ij}\) values in the each row add up to \(\lambda\), while remaining close to the original \(y_{ij}\) values.

globalMean <- mean(Y)
rowMean <- rowMeans(Y)
vecModel <- as.numeric(apply(Y, 2, function(x) x + (globalMean - rowMean)))
nCompNeg <- sum(vecModel < 0)
if (nCompNeg == 0) 
  Z2 <- matrix(vecModel, nrow = n, ncol = q)

However, these estimates can be negative. Back again to the full model which combines equality and inequality constraints, in general the estimates must be only numerically obtained. This estimation problem is a particular case of the optimization problem known as quadratic programing.

Proofs

Quadratic programming optimizes the objective \[ \frac{1}{2}\mathbf{b}^\textrm{T}\mathbf{D}\mathbf{b}−\mathbf{d}^\textrm{T}\mathbf{b} \quad \text{subject to }\,\mathbf{A}^\textrm{T} \ge \mathbf{b}_0 \]

However, we have

\[ \min_{\lambda_{ij}} \; \sum_{i=1}^{n} \sum_{j=1}^q \left( y_{ij} - \lambda_{ij} \right)^{2} = \min_{\Lambda}\left\| \mathbf{Y}-\mathbf{\Lambda} \right\|_F^2 \]

If we vectorize \(\mathbf{Y}\) and \(\mathbf{\Lambda}\)

\[ \mathbf{y} = \text{vec}(\mathbf{Y}) \quad\quad \boldsymbol{\lambda} = \text{vec}(\mathbf{\Lambda}) \]

we have

\[ \| \mathbf{Y} - \mathbf{\Lambda} \|_F^2 = \| \mathbf{y} - \boldsymbol{\lambda} \|_2^2 = (\mathbf{y} - \boldsymbol{\lambda})^\textrm{T} (\mathbf{y} - \boldsymbol{\lambda}) = \mathbf{y}^\textrm{T} \mathbf{y} - 2 \mathbf{y}^\textrm{T} \boldsymbol{\lambda} + \boldsymbol{\lambda}^\textrm{T} \boldsymbol{\lambda} \] \(\mathbf{y}^\textrm{T} \mathbf{y}\) is constant, so the optimization problem is equivalent to

\[ \min_{\boldsymbol{\lambda}} \; \boldsymbol{\lambda}^\textrm{T} \boldsymbol{\lambda} - 2 \mathbf{y}^\textrm{T} \boldsymbol{\lambda} = \min_{\boldsymbol{\lambda}} \; \frac{1}{2}\boldsymbol{\lambda}^\textrm{T}\mathbf{I}\boldsymbol{\lambda} - \mathbf{y}^\textrm{T} \boldsymbol{\lambda} \] where \(\mathbf{I}\) is the \(nq \times nq\) identity matrix.

if (nCompNeg > 0) {
  D <- diag(n*q)
  d <- as.numeric(t(Y))
  A <- matrix(0, nrow = n * q, ncol = (n * q + n - 1))
  A[1:q, 1:(n - 1)] <- 1
  for (i in 1:(n * q)) {
      A[i, (n - 1 + i)] <- 1
  }
  for (i in 2:n) {
      A[(((i - 1) * q) + 1):(i * q), (i - 1)] <- -1
  }
  sQP <- quadprog::solve.QP(D, d, A, meq = (n - 1))
  Z2 <- matrix(sQP$solution, nrow = n, ncol = q, byrow = T)
}

colnames(Z2) <- colnames(Y)
rownames(Z2) <- rownames(Y)

Schematic representation of the restriction matrix

The model of constant size projects the original observations on the surface of the \((q-1)\)-simplex. The projected data is \(\mathbf{Z}_{(2)}\).

\(V_{ls(2)}\) is the variability of the projected points:

\[ V_{ls(2)} = tr\left(\mathbf{\Sigma_{Z(2)}}\right) \tag{8}\]

sigma.z2 <- cov(Z2)
sigma.z2 <- sigma.z2 * (n - 1)
ss.t.z2 <- sum(diag(sigma.z2)) 

The ratio \(V_{ls(2)}/V_t =\) 0.1347664 represents the fraction of total variability explained by a constant size model.

Calculation of \(V_{ls(2)}\) (II)

OLS with restrictions via quadratic programming may be computationally expensive for large datasets. Hence, we propose an alternative optimization strategy.

Let’s consider a simplex of size \(\lambda\) (see Equation 6). We first make an orthogonal projection of the samples onto the hyperplane \(\{ \mathbf{x} \in \mathbb{R}^q \mid \mathbf{1}^\top\mathbf{x} = \lambda \}\). This is equivalent to row-centering \(\mathbf{Y}\) (orthogonal projection of row vectors on the hyperplane \(\{ \mathbf{x} \in \mathbb{R}^q \mid \mathbf{1}^\top\mathbf{x} = 0 \}\)) followed by adding the simplex center \(\mathbf{c} = (\frac{\lambda}{q}, \ldots, \frac{\lambda}{q})\) to each row:

\[ \mathbf{Z}_{(2)} = \mathbf{YC} + \frac{\lambda}{q}\mathbf{J}_n \tag{9}\]

where \(\mathbf{C} = \mathbf{I} - \frac{1}{n} \mathbf{J}_n\) and \(\mathbf{J}_n = \mathbf{1} \mathbf{1}^\top\).

Then, we pull the projections that have fallen outside the simplex back inside it (we enforce the constraint \(z_{(2)j} \ge 0 \, \forall j = 1, \ldots, q\)), using Duchi’s algorithm for Euclidean projection of a vector onto the simplex.

Duchi’s algorithm

---

INPUT: A vector \(\mathbf{v} \in \mathbb{R}^n\) and a scalar \(z > 0\).

\(\quad\) Sort \(\mathbf{v}\) into \(\boldsymbol{\mu}\): \(\mu_1 \ge \mu_2 \ge \cdots \ge \mu_n\)

\(\quad\) Find \(\rho = \max \left\{ j \in [n] : \mu_j - \frac{1}{j}\left( \sum_{r=1}^j \mu_r - z \right) > 0 \right\}.\)

\(\quad\) Define \(\theta = \frac{1}{\rho} \left(\sum_{i=1}^{\rho} \mu_i - z \right)\)

OUTPUT: \(\mathbf{w}\) s.t. \(w_i = \max \{ v_i - \theta, \, 0 \}\)

---

There are alternative formulations and algorithms for the same purpose (e.g. Wang & Carreira, Chen & Ye or Condat).

Framed this way, the optimization problem becomes:

\[ \min_{\lambda > 0} \left\| \mathbf{Y} - \mathbf{Z}_{(2)} \right\|_F^2 \] Note that for a given \(q\), \(\mathbf{Z}_{(2)}\) only depends on \(\lambda\) (see Equation 9). Hence, this is no longer an optimization in many variables but in just one (the size of the simplex, \(\lambda\)).

projectSimplex <- function(data, size = 1){
  
  # Enforce 1^Ty = size
  center <- rep(size/q, q)
  C <- diag(q) - matrix(1/q, q, q)
  Z <- data%*%C + center  

  # Enforce z(lambda)j >= 0 for points out of the simplex
  duchi <- function(v, size = 1) {
    s <- sort(v, decreasing = TRUE)
    cssv <- cumsum(s)
    rho <- max(which((s - (cssv - size) / seq_along(cssv)) > 0))
    theta <- (sum(s[1:rho]) - size)/rho
    w <- pmax(v - theta, 0)
    return(w)
  }
  
  out <- rowSums(Z < 0) == 1 # out of simplex?
  
  if (any(out)) {
    outSimplex <- Z[out, , drop = FALSE]
    outBackIn <- t(apply(outSimplex, 1, duchi, size))
    Z[out, ] <- outBackIn
  }
  
  return(Z) # rowSums(Z) = 1; all(Z >= 0) = TRUE
}

ssdifSimplex <- function(data, size = 1) {
    proj_data <- projectSimplex(data, size)
    return( sum((data - proj_data)^2) )
}

resOpt <- optimize(f = ssdifSimplex, 
                   interval = c(0, max(rowSums(Y))), # search space
                   data = Y
                   )

Z2b <- projectSimplex(Y, size = resOpt$minimum)

sigma.z2b <- cov(Z2b)
sigma.z2b <- sigma.z2b * (n - 1)
ss.t.z2b <- sum(diag(sigma.z2b))

Finally, we obtain \(\mathbf{Z}_{(2)}\) for the value of \(\lambda\) that minimizes our new objective and calculate \(V_{ls(2)}\) using Equation 8. Then we compute \(V_{ls(2)}/V_t =\) 0.1347664. Note that we obtain the same result as with the quadratic programming optimization strategy.

The search space for \(\lambda\) is highlighted.

Covariates

See the section on Covariates above. Essentially, now we need to project on the simplex the fitted values from the linear model and compute the corresponding sums of squares.

Y.lm.vls2 <- projectSimplex(Y.lm, size = resOpt$minimum) # alternatively using QP
ss.b.z2 <- sum(diag(cov(Y.lm.vls2))) * (n - 1) # ss.b already computed

The ratios of interest are \(\frac{SS_b}{SSt} =\) 0.275994 (the fraction of total variability explained by the covariate, already computed above) and \(\frac{SS_b(\mathbf{Z}_{(2)})}{SS_b} =\) 0.236794 (the fraction of the variability explained by the covariate that is explained by a constant size model).

Proposed updates

1. New R package
2. Evaluate and optimize performance
3. Outlier samples (bootstrap?)
4. Incorporate inference