(This response picks up where @AVB, who has provided useful comments, left off by suggesting we need to figure out which differences $X_i - X_j$ ought to be included among the independent variables.)
The big question here is what is an effective method to identify the model. Later we can worry about faster methods. (But regression is so fast that you could process dozens of variables for millions of records in a matter of seconds.)
To make sure I'm not going astray, and to illustrate the procedure, I simulated a dataset like yours, only a little simpler. It consists of 60 independent draws from a common multivariate normal distribution with five unit-variance variables $Z_1, Z_2, Z_3, Z_4,$ and $Y$. The first two variables are independent of the second two and have correlation coefficient 0.9. The second two variables have correlation coefficient -0.9. The correlations between $Z_i$ and $Y$ are 0.5, 0.5, 0.5, and -0.5. Then--this changes nothing essential but makes the data a little more interesting--I rescaled the variables, thus: $X_1 = Z_1, X_2 = 2 Z_2, X_3 = 3 Z_3, X_4 = 4 Z_4$.
Let's begin by establishing that this simulation emulates the stated problem. Here is a scatterplot matrix.
The full regression of $Y$ against the $X_i$ is highly significant ($F(4, 55) = 15.28,\ p < 0.0001$) but all four t-values equal 1.24 ($p = 0.222$), which is not significant at all. The estimated coefficients are 0.26, 0.13, 0.088, and -0.066 (rounded to two sig figs).
Here is my proposal: systematically combine variables in pairs (six pairs in this case, 36 pairs for nine variables), one pair at a time. Regress a pair along with all remaining variables, seeking highly significant results for the pairs.
What is a "pair"? It is the linear combination suggested by the estimated coefficients. In this case, they are
$$\eqalign{
X_{12} =& X_1 / 0.26 &+ X_2 / 0.13 \cr
X_{13} =& X_1 / 0.26 &+ X_3 / 0.088 \cr
X_{14} =& X_1 / 0.26 &- X_4 / 0.066 \cr
X_{23} =& X_2 / 0.13 &+ X_3 / 0.088 \cr
X_{24} =& X_2 / 0.13 &- X_4 / 0.066 \cr
X_{34} =& X_3 / 0.088 &- X_4 / 0.066 \text{.}
}$$
In general, with $\hat{\beta}_i$ representing the estimated coefficient of $X_i$ in this full regression, the pairs are defined by
$$X_{ij} = X_i / \hat{\beta}_i + X_j / \hat{\beta}_j\text{.}$$
This is so systematic that it's straightforward to script.
The "identification regressions" are the model
$$Y \sim X_{12} + X_3 + X_4$$
along with the five additional permutations thereof, one for each pair.
You are looking for results where $X_{ij}$ becomes significant: ignore the significance of the remaining $X_k$. To see what's going on, I will list the results of all six identification regressions for the simulation. As a shorthand, I list the variables followed by a vector of their t-values only:
$$\eqalign{
X_{12}, X_3, X_4:&\ (5.50, 1.24, -1.24) \cr
X_{13}, X_2, X_4:&\ (1.36, 4.94, -1.13) \cr
X_{14}, X_2, X_3:&\ (1.31, 5.16, 1.17) \cr
X_{23}, X_1, X_4:&\ (1.64, 3.10, -1.09) \cr
X_{24}, X_1, X_3:&\ (1.50, 4.15, 1.07) \cr
X_{34}, X_1, X_2:&\ (5.56, 1.25, 1.25)
}
$$
As you can see from the first component of each vector (the t-value for the pair), precisely two disjoint pairs exhibit significant t-statistics: $X_{12}$, with $t = 5.50\ (p \lt 0.001)$, and $X_{34}$, with $t = 5.56\ (p \lt 0.001)$. The model thus identified is
$$Y \sim X_{12} + X_{34}\text{.}$$
(In general, we would also include--provisionally--any remaining $X_i$ not participating in any of the pairs. There aren't any in this case.)
The regression results are
$$\eqalign{
\hat{\beta_{12}} &= 0.027\ (t = 5.54,\ p \lt 0.001) \cr
\hat{\beta_{34}} &= 0.0055\ (t = 5.58,\ p \lt 0.001), \cr
F(2, 57) &= 30.92\ (p \lt 0.0001).
}$$
Translating back to the original $X_i$, the model is
$$\eqalign{
Y &= 0.027(X_1 / 0.26 + X_2 / 0.13) + 0.0055(X_3 / 0.088 - X_4 / 0.066) \cr
&= 0.103 X_1 + 0.206 X_2 + 0.0629 X_3 - 0.0839 X_4 \cr
&= 0.103 (Z_1 + Z_2) + 0.021 (Z_3 - Z_4) \text{.}
}$$
(The last line shows how this all relates to form of the original question.) That's exactly the form used in the simulation: $Z_1$ and $Z_2$ enter with the same coefficient and $Z_3$ and $Z_4$ enter with opposite coefficients. This method got the right answer.
I want to share a cool observation in this regard. First, here's the scatterplot matrix for the model.
Notice how $X_{12}$ and $X_{34}$ look uncorrelated. Furthermore, $Y$ is only weakly correlated with these variables. Doesn't look like much of a relationship, does it? Now consider an alternative set of pairs, $X_{13}$ and $X_{24}$. The regression of $Y$ on these is still highly significant ($F(2, 57) = 16.61\ (p \lt 0.0001).$ Moreover, the coefficient of $X_{24}$ is significant ($t = 2.39,\ p = 0.020$) even though that of $X_{13}$ is not ($t = 0.24,\ p = 0.812$). But look at the scatterplot matrix!
Clearly $X_{13}$ and $X_{24}$ are strongly correlated. But, even though this is the wrong model, $Y$ is also visibly correlated with these two variables, much more so than in the preceding scatterplot matrix!
The lesson here is that mere bivariate plots can be deceiving in a multiple regression setting: to analyze the relationship between any candidate independent variable (such as $X_{12}$) and the independent variable ($Y$), we must make sure to "factor out" all other independent variables. (This is done by regressing $Y$ on all other independent variables and, separately, regressing $X_{12}$ on all the others. Then one looks at a scatterplot of the residuals of the first regression against the residuals of the second regression. It's a theorem that the slope in this bivariate regression equals the coefficient of $X_{12}$ in the full multivariate regression of $Y$ against all the variables.)
This insight shows why we might want to systematically perform the "identification regressions" I have proposed, rather than using graphical methods or attempting to combine many of the pairs in one model. Each identification regression assesses the strength of the contribution of a proposed linear combination of variables (a "pair") in the context of all the remaining independent variables.
Note that although correlated variables were involved, correlation is not an essential feature of the problem or of the solution. Even where you don't expect the original variables $X_i$ to be strongly correlated, you could expect a model to have (unknown) linear constraints among the variables. That is the important issue to cope with. The presence of correlation only means that it can be problematic to identify such pairs solely by inspecting the original regression results.
Following the procedure I have proposed does not guarantee you will find a unique solution. It's conceivable, for instance, that you will find so many highly significant pairs that they are linearly dependent, forcing you to select among them by some other criterion. Nevertheless, the results you get ought to limit the sets of pairs you need to examine; they can be obtained with a straightforward procedure without intervention; and--if this simulation is any guide--they have a good chance of producing effective results.
