# Sampling Normal variables with linear constraints and given variances - Fraser (1951)

I am trying to generate $$\mathbf{X} = (X_1, \ldots, X_d)$$ such that $$X_i \sim N(0, v_i)$$ for $$i = 1, \ldots, d$$ satisfying the linear constraint $$\sum_{j = 1}^dX_j = 0$$.

The work of Fraser (1951) proposes a solution to this problem. Although I was able to understand it, I was not able to implement it completely. My steps are below:

Steps:

1. Rewrite the linear constraint $$\sum_{j = 1}^d X_j = 0$$ as $$\sum_{j = 1}^d b_{dj} X_j = 0$$ with $$||\mathbf{b}_d||_2 = 1$$ which yields: $$\mathbf{b}_d = [1/\sqrt{d}, \ldots, 1/\sqrt{d}].$$

2. Add d-1 rows to complete a matrix $$\mathbf{B}_{d \times d}$$ satisfying orthogonality conditions: $$<\mathbf{b}_k, \mathbf{b}_l> = \begin{cases} 0 \quad \mbox{if k \neq l} \\ 1 \quad \mbox{if k = l}. \end{cases}$$

Based on Gramm-Schmidt process, the $$i^{th}$$ row of $$\mathbf{B}_{d \times d}$$ - denoted as $$\mathbf{b}_i$$ - is defined as $$b_{ij} \begin{cases} 0 \quad &\mbox{if j < i} \\ \left(\displaystyle\frac{d - i}{d - i + 1}\right)^{1/2} \quad& \mbox{if i = j} \\ - \displaystyle\frac{1}{[(d - i + 1)(d - i)]^{1/2}} \quad& \mbox{if j > i}. \end{cases}$$

Double-checking in R:

d <- 4

B <- matrix(NA, nrow = d, ncol = d)
B[d, ] <- rep((1/sqrt(d)), d)

for (i in 1:(d-1)){
for (j in 1:d){
if (j < i)
B[i, j] <- 0
if (i == j)
B[i, j] <- ((d - i)/(d - i + 1))^(1/2)
if (j > i)
B[i, j] <- - 1/(((d - i + 1)*(d-i))^(1/2))
}
}

# Normalized
k <- 1
l <- 1
B[k, ] %*% B[l, ]
#     [,1]
#[1,]    1

# Orthogonal
k <- 1
l <- 2
B[k, ] %*% B[l, ]
#              [,1]
#[1,] -3.469447e-18


The matrix $$\mathbf{B}$$ can be considered as an orthogonal rotation of $$n$$-space. Consider coordinates $$\mathbf{Y} = (Y_1, \ldots, Y_n)$$ with respect to the new axes:

$$Y_i = \sum_{j = 1}^d b_{ij} X_j.$$

In a matrix notation, $$\mathbf{Y} = \mathbf{B}\mathbf{X}$$. A set of constrained variables $$(X_1, \ldots, X_n)$$ is transformed to a set of unconstrainted variables $$(Y_1, \ldots, Y_{n-1}, Y_n = 0)$$.

Double-checking in R:

X <- c(1, 2, 3, -6)
Y <- B%*%X
Y
#         [,1]
#[1,] 1.154701
#[2,] 2.857738
#[3,] 6.363961
#[4,] 0.000000

1. Establish variances of the variables $$\mathbf{Y}$$ based on the distribution of $$\mathbf{X}$$. Given that $$\mathbf{X} \sim N_d(\mathbf{0}, \mathbf{V})$$ with $$\mathbf{V} = diag(v_1, \ldots, v_d)$$, it follows that

$$\mathbf{Y} = \mathbf{B}_{(d-1)\times(d-1)}\mathbf{X} \sim N_{d-1}(\mathbf{0}, \boldsymbol{\tau}),$$ where $$\boldsymbol{\tau}_{(d-1)\times(d-1)} = \mathbf{B}_{(d-1)\times d}\mathbf{V}_{d\times d)}\mathbf{B}_{(d-1)\times d}^T$$ with $$\mathbf{A}_{l\times k}$$ is a sub matrix with the first $$l$$ rows and $$k$$ columns of $$\mathbf{A}$$.

In our running example assuming that $$v_i = 1 \forall i$$ for simplicity,

V <- diag(1, nrow = d, ncol = d)
tau <- B[1:(d-1), ]%*%V%*%t(B[1:(d-1), ])
> tau
#              [,1]         [,2]          [,3]
#[1,]  1.000000e+00 2.775558e-17 -2.775558e-17
#[2,]  2.775558e-17 1.000000e+00  0.000000e+00
#[3,] -2.775558e-17 0.000000e+00  1.000000e+00

1. Generate $$\mathbf{Y}$$ as described in step 4. Then, calculate $$\mathbf{X} = \mathbf{B}^{-1}\mathbf{Y} = \mathbf{B}^{T}\mathbf{Y}$$ because $$\mathbf{B}$$ is an orthonormal matrix.
set.seed(1234)
aux <- matrix(rnorm((d-1)*100000,0,1),100000,d-1)
y <- cbind(aux%*%chol(tau), 0)

x <- y %*%B


However, this procedure does not produce the desirable variances:

mean(rowSums(x))
# [1] -9.463654e-19
apply(x, 2, var)
# [1] 0.7491879 0.7423549 0.7546185 0.7540128


I think that my error is probably at step 3 since the equality $$\tau_{ii} = v_i$$ for $$i = 1, \ldots, d$$ always will hold because $$\mathbf{B}$$ is an orthonormal matrix. However, my intuition is that $$\tau_{ii}$$ should be higher than $$v_i$$. I could not grasp very well the discussion in Fraser's section 4 summarized as steps 5.3 and 5.4 in section 5.

Any thoughts?

• Your mistake is in $\tau$, you cannot use the normality assumption to state that $\boldsymbol{\tau} = \boldsymbol{B}\boldsymbol{V}\boldsymbol{B}^{\prime}$. I will write a solution Commented Dec 28, 2021 at 3:58

Suppose $$(Y_1, \cdots, Y_{n-1}) \sim N_{n-1}\left(\boldsymbol{0}_{n-1}, \boldsymbol{T}\right)$$. For simplicity, let $$v_i=1$$ for all $$i$$ and let $$b_i^{\ast}$$ denote the $$i$$th column of $$\boldsymbol{B}$$ excluding observation $$b_{ni}$$. That reference states that the elements of $$\boldsymbol{T} = (\tau_{rs})$$ are found by solving the set of $$n$$ equations $$\begin{eqnarray*} \sum_{r=1}^{n-1}\sum_{s=1}^{n-1} b_{ri}\tau_{rs}b_{si} = 1 \quad \mbox{for} \quad i=1,\cdots,n. \end{eqnarray*}$$ Note that this can be expressed as the equivalent quadratic form $$\left(b_i^{\ast}\right)^{\prime}\boldsymbol{T}b_i^{\ast}=1.$$

Now we wish to solve for the elements of $$\boldsymbol{T}$$. Since the quadratic form is a scalar, it is equal to its trace. Using the fact that $$\mbox{tr}(AB)=vec^{\prime}(A^{\prime})vec(B)$$, $$\mbox{tr}(AB)=\mbox{tr}(BA)$$, and for a symmetric $$n \times n$$ matrix, $$A$$, that $$D_n vech(A)=vec(A)$$, with $$vec(\cdot)$$ denoting the vectorization of a matrix, $$vech(\cdot)$$ denoting the half-vectorization of a matrix and $$D_n$$ denoting the duplication matrix of order $$n$$. Since $$vech(\boldsymbol{T})$$ contains all of the unique elements of $$\boldsymbol{T}$$, we can create a linear system of equations using the equivalent formulation: $$\begin{eqnarray*} \left(b_i^{\ast}\right)^{\prime}\boldsymbol{T}b_i^{\ast} &=& \mbox{tr} \left(\left(b_i^{\ast}\right)^{\prime}\boldsymbol{T}b_i^{\ast}\right) \\ &=& \mbox{tr} \left(b_i^{\ast}\left(b_i^{\ast}\right)^{\prime}\boldsymbol{T}\right) \\ &=& vec^{\prime}\left(b_i^{\ast}\left(b_i^{\ast}\right)^{\prime}\right)vec\left(\boldsymbol{T}\right) \\ &=& vec^{\prime}\left(b_i^{\ast}\left(b_i^{\ast}\right)^{\prime}\right) \boldsymbol{D}_{n-1} vech\left(\boldsymbol{T}\right). \end{eqnarray*}$$

Next let $$\begin{eqnarray*} \boldsymbol{W} = \begin{pmatrix} vec^{\prime}\left(b_1^{\ast}\left(b_1^{\ast}\right)^{\prime}\right) \boldsymbol{D}_{n-1} \\ \vdots \\ vec^{\prime}\left(b_n^{\ast}\left(b_n^{\ast}\right)^{\prime}\right) \boldsymbol{D}_{n-1}. \end{pmatrix} \end{eqnarray*}$$ Hence the original $$n$$ set of equations can be written as $$\begin{eqnarray*} \boldsymbol{W} vech\left(\boldsymbol{T}\right) = \boldsymbol{1}_n, \end{eqnarray*}$$ and for any generalized inverse of $$\boldsymbol{W}^{\prime}\boldsymbol{W}$$, we have $$\begin{eqnarray*} vech\left(\boldsymbol{T}\right) = \left(\boldsymbol{W}^{\prime}\boldsymbol{W}\right)^{-} \boldsymbol{W}^{\prime}\boldsymbol{1}_n \end{eqnarray*}$$.

In fact, for this problem, one may reduce the above to find that $$\boldsymbol{T} = \frac{n}{n-1} \boldsymbol{I}_{n-1}$$.

Using the $${\tt matrixcalc}$$ package in $${\tt R}$$ to obtain the duplication matrix, the correct code should be

d <- 4

B <- matrix(NA, nrow = d, ncol = d)
B[d, ] <- rep((1/sqrt(d)), d)

for (j in 1:(d-1)){
for (i in 1:d){
if (i < j)
B[j, i] <- 0
if (i == j)
B[j, i] <- ((d - j)/(d - j + 1))^(1/2)
if (i > j)
B[j, i] <- - 1/(((d-j+1)*(d-j))^(1/2))
}
}

set.seed(1234)
library(matrixcalc)
library(Matrix)
library(MASS)
D = duplication.matrix(d-1)
W = matrix(0,d,d*(d-1)/2)
for (i in 1:d){
W[i,] = t(as.vector(B[-d,i]%*%t(B[-d,i])))%*%D
}
vech.tau = as.vector(ginv(t(W)%*%W)%*%t(W)%*%rep(1,d))
tau = matrix(0,d-1,d-1)
l=1
for(i in 1:(d-1)){
for (j in i:(d-1)){
tau[i,j] = vech.tau[l]
l = l+1
}
}
tau = forceSymmetric(tau)
aux <- matrix(rnorm((d-1)*100000,0,1),100000,d-1)
sqrt_tau <- chol(tau)
y <- cbind(aux%*%chol(tau), 0)

x <- y %*%B
mean(rowSums(x))
# [1] -6.543378e-19
apply(x, 2, var)
# [1] 0.9989172 0.9898065 1.0061580 1.0053505
$$$$
`