Linear regression where attribute vectors can be reordered 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.
 A: 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)$$
where


*

*$\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.
#
set.seed(17)
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),
                           method="BFGS"))
#
# The fitted coefficients should be compared to `beta` *up to a permutation*!
#
model
sqrt(model$value / (n-1)) # (Under)estimates the SD of `eps` $
sd(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)
    model$value

