Multiple linear regression for multi-dimensional input and output? Assume that I have $N$ points $x_i,i=1,...,N$ in some $A&gt;1$-dimensional  space $\mathbb{R}^A$ with pointwise evaluations of some function $f:\mathbb{R}^A \rightarrow \mathbb{R}^B$, i.e. $f(x_i),i=1,...,N$ where $f(x_i) \in \mathbb{R}^B$.
It is my goal to find a multiple linear regression between $x_i$ and $f(x_i)$. Now sklearn has a function (sklearn.linear_model.LinearRegression) for a multiple linear regression for functions of the type $f:\mathbb{R}^A \rightarrow \mathbb{R} $, but my output is $B$-dimensional. I assume that I could make independent multiple linear regressions for each output dimension and then combine the results, but there must be an easier way of achieving this.
Do you know of a more efficient way?
You are asking about multioutput regression. The class you talked about sklearn.linear_model.LinearRegression supports this out of the box.
import numpy as np
from sklearn.linear_model import LinearRegression

# features
A = 10
# number of values to predict
B = 15
# number of rows in dataset
m = 100

x = np.ones((m, A))
y = np.ones((m, B))

model = LinearRegression()
model.fit(x, y)

sklearn.linear_model.LinearRegression actually just creates B models. However it optimises calculations using vectorisation.
It actually is exactly the same as a fully connected layer in a neural network which has no activation function.
You can read more about it here:
https://machinelearningmastery.com/multi-output-regression-models-with-python/