I'm attempting to do a constrained multi-variate regression of the form


Where Y is a Nx24 matrix, X is a Nx1 vector, B is a 1x24 vector and E is a Nx24 matrix

I can solve the constrained multiple regression problem using say solve.QP or lsei. However it appears all these solvers expect Y to be a vector. For example the unconstrained case can be solved using limSolve::Solve


A <- as.matrix(X)
B <- Y
s <- Solve(A,B)

But the method used to solve the generalized form

lsei(A = A, B = B, fulloutput = TRUE, verbose = FALSE)

complains that A and B are incompatible, and the documentation specifically states that B is a vector. Similarly solve.QP tests the length of dvec which only works (correctly at least) for a vector.

Is there's a generalised solver in R that can solve this problem? Or am I using the standard solvers incorrectly.


If you are going to use one of those standardized optimization solvers "off the shelf", you will need to "vec" your matrices by stacking the columns such that your matrix of parameters and LHS becomes a single vector.

However, an easier and less error-prone way is to use a modeling tool which will accept your problem in original (in this case, matrix) form, and do the necessary conversions under the hood to reformulate the problem in a way which the optimization solver it calls can handle, and then transform the solver's results back into your original form. Given that you are using R, I recommend you consider using CVXR https://cran.r-project.org/web/packages/CVXR/index.html , which should easily handle your problem, which is right in its wheelhouse, and give you easy growth opportunity to handle many variants of this problem as you find the need. If you read the beginning of https://cvxr.rbind.io/post/examples/cvxr_gentle-intro/ , you ought to be in business very quickly.

  • 1
    $\begingroup$ Thanks, I'll take a look at CVXR - it seems to be what I'm after. $\endgroup$ – David Waterworth May 24 '18 at 23:24
  • 1
    $\begingroup$ I actually ended up stacking, but CVXR looks excellent, thanks for pointing me towards it. I'll most definitely be using it in the future. $\endgroup$ – David Waterworth May 26 '18 at 2:26
  • 1
    $\begingroup$ If you do use CVXR and have LHS - RHS consisting of a matrix, then you would want to use Frobenius norm as the objective (no need to square, and better not to): norm(LHS - RHS, "F"), which winds up being the same as the 2-norm of the vec'd (stacked version) of LHS - RHS. $\endgroup$ – Mark L. Stone May 26 '18 at 2:32

I think most implementations try to fit a matrix to a vector. $\vec{y}$ would contain your observations and $X = (\vec{x_1}, \vec{x_2}, \dots)$ vectors of input features. So, it is a matrix... then you have the weight matrix $W$ and finally your bias vector $\vec{b}$.

$$ W\cdot X + \vec{b} = \vec{y}$$

Then you fit weights and biases by minimize a certain loss function. E.g. mean squared error. But in your case it must be a distance between matrices. I suppose, most standard implementation don't handle that. Just check the documentation strings of your functions. Usually, input and output requirements are specified.

You should either tweak your problem, e.g. by predicting each row in Y separately or implement your own function, which should not be too difficult. ...or search for an implementation that generalizes to your case.


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