What is the proper way to do vector based linear regression in R I want to do linear regression between vector inputs and vector output. That is each $y$ is a vector with $M$ components, and each $x$ is a vector with $N$ components and the answer should look like $y \sim Ax + b$ where $A$ is an $M \times N$ matrix and $b$ is a vector with $M$ components.
I have a very clear understanding of the concept and what I want R to do, but it is the proper syntax I am lacking.
Trying to google around to find this has been quite difficult because terms like multivariable seem to always point me to answers of the form 
$$y \sim x_1 + x_2 + x_3 + \dots + x_n$$
where there are multiple input sources (or rather, a multidimensional input), but never with multidimensional outputs.
If I just feed in matrices for $y$ and $x$ that MIGHT give what I want, but it might also just treat each $y$ component as directly related to each $x$ component and give an answer based on that ($M = N$ for the important instance I have). So I have to be sure that I am doing it correctly.
What is the correct means for using R to do linear regression of the sort
$$y \sim A x + b $$
where the solution $A$ is an $M \times N$ matrix, and $b$ is a vector of length $M$, and each datum $x$ is a vector of length $N$ and each corresponding datum $y$ is a vector of length $M$?
 A: Set up the problem. We have 4 "data" items
guys=c("fred","bob","andrew","joe");
## [1] "fred"   "bob"    "andrew" "joe"   

Each row $x$ of $X$ is a vector with $N=3$ components
# By default, matrices are filled by column
X=matrix(c(1,0,1,0,  1,0,0,1,  1,1,1,0), ncol=3,
  dimnames=list(guys,c("tea","coffee","yoga")))
##        tea coffee yoga
## fred     1      1    1
## bob      0      0    1
## andrew   1      0    1
## joe      0      1    0

Each row $y$ of $Y$ is a vector with $M=2$ components
Y=matrix(c(185,209,138,167,  36, 32, 30, 31), ncol=2,
  dimnames=list(guys, c("weight","girth")))
##        weight girth
## fred      185    36
## bob       209    32
## andrew    138    30
## joe       167    31

This is how we can fit a linear model $Y=X\times A$ with no intercept
fit = lm(Y~X+0)
## ...
## Coefficients:
##          weight   girth
## Xtea     -111.00   -10.33
## Xcoffee   127.00    22.67
## Xyoga     209.00    32.00

It is the same as calling lm.fit(X,Y)
fit_ = lm.fit(X,Y)
sum(abs(coef(fit)-coef(fit_)))
## 0

except that you get a decent "summary()" etc.
summary(fit)
## ...
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)
## Xtea      -10.33      18.63  -0.555    0.678
## Xcoffee    22.67      11.79   1.923    0.305
## Xyoga      32.00      14.43   2.217    0.270

Check that the coefficient matrix has the intended meaning of $A$ in $Y=X\times A$
sum(abs(
  X %*% coef(fit) - Y
))
## [1] 145
sum(abs(residuals(fit)))
## [1] 145

Here $A$ is an $N \times M$ matrix, while you asked for $M \times N$. If you really wanted $Y=A \times X$ instead of $Y=X \times A$, then take the transpose of everything. Each "guy" will then be a column rather than a row. But usually in R we think of observations being stored in rows, that's why I did it this way.
Here is the version with an intercept
fit_int = lm(Y~X+1)     # same as default: lm(Y~X)
##              weight  girth
## (Intercept)  120      25
## Xtea         -71      -2
## Xcoffee       47       6
## Xyoga         89       7

The intercept is part of the "coef()" matrix, so instead of $Y=X\times A+b$ we
have to modify $X$, it's like $Y=[1,X]\times A$:
sum(abs(
  cbind(1,X) %*% coef(fit_int) - Y
))
## [1] 2.842171e-14
sum(abs(residuals(fit_int)))
## [1] 0   # same as above, more or less

In my setup, after introducing an intercept term, the number of unknowns now matches the number of constraints and the equation has an exact solution as you can see.
I added row and column names for clarity, unfortunately R does not check that these match when doing matrix arithmetic or when creating a formula. So be careful you don't get rows and columns mixed up, or columns out of order. The following produces no error message for example:
t(cbind(X,1)) %*% coef(fit_int)  # wrong!

