Is there an existing method to solve this kind of regression problem?

Given a set of $n$ training data $\{(x_1,y_1)\ldots(x_n,y_n)\}$, find $\beta$ and a set of functions $\{f_1,f_2,\ldots,f_n\}$ to minimize $\sum_{i}\left(y_i-f_i(x_i)^T\beta\right)^2$ where each $x_i$ is a length $m$ vector, $\beta$ is a length $m$ weight vector, and $f_i$ is a function that permutes the elements of $x_i$. For example, if $x_i=(a,b,c)^T$, then $f_i(x_i)$ is one of the six possible permutations $(a,b,c)^T,(a,c,b)^T,(b,a,c)^T,(b,c,a)^T,(c,a,b)^T,(c,b,a)^T$.

Note that without the permutation function $f_i$, the objective function is just the residual sum of squares, $\sum_{i}\left(y_i-x_i^T\beta\right)^2$, and this becomes a typical linear regression problem.

This seems like a hard problem. It would be too computationally expensive to solve it by brute force, e.g. perform a normal linear regression on every possible reordering of every possible data point.

The only saving grace is that I know a lot of the reorderings will be the same, e.g. the first 200 data points will have the same ordering, the second 200 data points will have the same ordering but possibly different from the first, etc.

This kind of problem comes up if there isn't any apparent way to order the attributes in your data. For example, each data point may represent the incomes of a group of 3 people (i.e. $x_i = (income_1,income_2,income_3)$), and each data point comes from a different group. If you don't know a reasonable property by which to associate members across groups, you can't order the incomes in a meaningful way.

  • 1
    $\begingroup$ What is the motivation for such a question? $\endgroup$ – Glen_b Sep 5 '13 at 8:22
  • $\begingroup$ Are you summing over all reorderings? Even computing the sum would be hard when n is non-trivial. But if it's over all reorderings I guess it could be possible to find some fast method for computing the sum. Similar to the sum over all subsets method in conditional logistic regression, described e.g. here: hydra.usc.edu/pm518b/materials%202012/… $\endgroup$ – Simen Gaure Sep 5 '13 at 9:18
  • 1
    $\begingroup$ I can't even understand the question because the notation is self-contradictory: expressions like "$a\in x_i$ imply $x_i$ is a set, not a vector; and $\beta f_i(x_i)$ is not defined because both $\beta$ and $f_i(x_i)$ are suppose to be elements of $\mathbb{R}^n$, which does not have a multiplication except for $n=1,2,4,8$. Do you think you could express your problem in words so that we don't have to figure out unconventional uses of mathematical expressions? $\endgroup$ – whuber Sep 5 '13 at 16:21
  • 1
    $\begingroup$ Thanks, it's clearer. I still wonder about one thing, though: in your example you use $3$-vectors $(a,b,c)$, whence $n=3$, but you list $6$ of the $f_i$, whence $n=6$. How is this contradiction to be resolved? I am also confused by what you are summing over: $i$ is the index, but now (with two different values of $n$) the term "$f_i(x_i)$" makes no sense. $\endgroup$ – whuber Sep 5 '13 at 20:36
  • 1
    $\begingroup$ @whuber: I don't think m will be larger than 7. $\endgroup$ – Joseph Sep 6 '13 at 4:27

For small numbers of variables, this problem is easier than it looks: for each group of data you can systematically try all $m!$ permutations of the coefficients and retain the permutation that minimizes the sum of residual squares.

To describe the algorithm in detail, let the data consist of $g$ pairs $X^{(j)}, Y^{(j)}$ where each $X^{(j)}$ is an $n_j$ by $m$ matrix whose columns are the independent values and each $Y^{(j)}$ is an $n_j$ column vector of dependent values. The model is

$$E[Y^{(j)}] = X^{(j)}\sigma_j(\beta)$$


  • $\beta$ is an $m$-vector of coefficients
  • $\sigma_j$ are permutations of the coefficients of $\beta$
  • $1 \le j \le g$.

Let $\mathfrak{S}_m$ denote the group of all permutations of the positions of an $m$-vector. To estimate the parameters $\beta$ and $\{\sigma_j\}$, define

$$f_j(\beta) = \min_{\sigma \in \mathfrak{S}_m} ||Y^{(j)} - X^{(j)}\sigma(\beta)||^2.$$

Given $\beta$, this amounts to finding a permutation of $\beta$ which minimizes the residuals in group $j$. One way to find this permutation is to compute $||Y^{(j)} - X^{(j)}\sigma(\beta)||^2$ for each of the $m!$ permutations $\sigma$. This brute force approach is reasonable for small $m$.

The estimate of $\beta$ minimizes the objective $F(\beta) = \sum_{j=1}^g f_j(\beta)$. Because $F$ is smooth almost everywhere, continuous, and quadratic in a neighborhood of is minimum, we may apply a gradient-based numerical solver to find the optimal value $ \widehat\beta$. This can be expected to converge quickly. Thus, the total effort is proportional (roughly) to $g m! \bar n$ where $\bar n = (n_1+\cdots+n_g)/g$ is the average group size.

Having obtained an estimate of $\beta$--which involves estimating the $\sigma_j$ when evaluating the $f_j$--we may then proceed as usual to predict $Y$, compute residuals, perform all diagnostics, etc. A convenient way to implement this would be to permute the rows of the data $X^{(j)}$ once and for all according to the inverse of the estimated $\sigma_j$ and then perform a standard regression on the dataset obtained by combining all the $X^{(j)}, Y^{(j)}$ cases.

The usual least-squares (and maximum likelihood) theory should hold approximately provided the groups are large and not too great in number. (When there is a large number of groups, some of the $\sigma_j$ are likely to be incorrectly estimated, leading to a reduction in the sum of squared residuals: in effect, we should expect to overfit the data.)

As a moderately difficult example I synthesized $1979$ cases with $m=6$ variables (each generated independently from a standard Normal distribution) partitioned into $g=10$ groups averaging almost $200 = \bar n$ per group. I set $\beta=(1,-2,3,-4,5,-6)$ and added iid Normally distributed error $\varepsilon$ with standard deviation $\tau=12$. After $50$ evaluations of $F$, requiring $15$ seconds total, convergence was achieved at the (reasonable) estimate $\widehat\beta = (1.10, -2.15, 2.88, -4.50, 4.98, -6.41)$. The estimate of $\tau$ was $12.12$ (to be compared to the actual realized value of $12.15$). Evidently overfitting was not a problem in this example.

Here is the R code used for the example. It is written in a moderately general way to allow flexible experimentation and even application to real data.

# Specify the problem.
m <- 6         # Number of variables (including any constant)
n.groups <- 10                            # Number of groups
n.per.group <- rbinom(n.groups, 250, 4/5)  # Group sizes
n <- sum(n.per.group)                      # Number of cases
eps <- rnorm(n, 0, 2*m)                    # Error terms
suppressWarnings(beta <- 1:m * c(1,-1))    # Coefficients
f <- t(replicate(n.groups, sample.int(m))) # Permutations
# Create data (y, x, g): IV, DVs, and groupings.
pred <- function(b, x, group, s) {
  # Apply the parameters (b, s) to (x, group).
  unlist(sapply(unique(group), function(g) x[group==g, ] %*% b[s[g, ]]))
group <- unlist(sapply(1:n.groups, function(i) rep(i, n.per.group[i])))
x <- matrix(rnorm(m*n), n, m)
dimnames(x) <- list(NULL, paste("x", 1:m, sep="."))
y <- pred(beta, x, group, f) + eps
# Define the objective function.
permutations <- function(x) {
  # Returns all the permutations of `x`, one per row.
  if (length(x) <= 1) return(matrix(x, 1, length(x)))
  matrix(sapply(1:length(x), function(i) t(cbind(permutations(x[-i]), x[i]))),
         ncol=length(x), byrow=TRUE)
fit <- function(y, x, b, sigma=permutations(1:length(b))) {
  # Returns a permutation `sigma` minimizing |y - x.sigma(b)|^2 along
  # with the minimal value.
  delta <- x %*% apply(sigma, 1, function(s) b[s])
  costs <- apply(delta, 2, function(z) {w <- z-y; sum(w*w)})
  i <- which.min(costs)
  list(permutation=sigma[i,], value=costs[i])
cost <- function(y, x, b, group=1, sigma=permutations(1:length(b))) {
  # Returns the total cost and the array of permutations minimizing it
  # for IVs `x`, DV `y`, and coefficients `b`.  The grouping of `x`
  # is given by `group`.  The feasible permutations are rows of `sigma`.
  sapply(unique(group), function(g) {
    i <- group==g
    fit(y[i], x[i, ], b, sigma)
  # Returns one column per group.
objective <- function(b, y, x, g, sigma=permutations(1:length(b))) {
  # Returns the value for coefficients `b` and IVs `x`, DV `y`, groups `g`.
  sum(unlist(cost(y, x, b, g, sigma)["value", ]))
# Find *some* starting value (albeit a poor one).
b <- coef(lm(y ~ x - 1))
b <- beta
# Optimize.
sigma <- permutations(1:length(b))
system.time(model <- optim(b, function(b) objective(b, y, x, group, sigma),
# The fitted coefficients should be compared to `beta` *up to a permutation*!
sqrt(model$value / (n-1)) # (Under)estimates the SD of `eps` $
# Diagnostics and checks.
f.hat <- t(sapply(cost(y, x, model$par, group, sigma)["permutation", ], identity))
    y.hat <- pred(model$par, x, group, f.hat)
pairs(cbind(x, y, y.hat), pch=".")
# Additional checks.
cost(y, x, b, group, sigma)
objective(beta, y, x, group, sigma)
objective(b, y, x, group, sigma)
objective(model$par, y, x, group, sigma)
  • $\begingroup$ I like this solution! One concern I have is that given so many possible permutations, wouldn't there be a lot of local minima to get stuck at? $\endgroup$ – Joseph Sep 6 '13 at 19:14
  • $\begingroup$ There could be if there exist any permutations that don't appreciably affect the fit in a group. This would tend to happen when the fit is awful in some group. If in fact there is some significant linear relationship for at least one permutation in every group, you shouldn't have any problems. Note that there will always be up to $m!$ equivalent local minima, corresponding to all the permutations of $\widehat\beta$. It doesn't matter which of these is found. $\endgroup$ – whuber Sep 6 '13 at 21:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.