I have a linear system where I am given 2 matrices, $A$ and $B$, and 2 vectors, $v$ and $c$, and I need to solve for the vector $x$. $A$ is $n\times n$, $B$ is $n \times n \times n$, and the vectors $v$, $c$, and $x$, are all $n\times 1$. I need to solve for $x$, which is complicated as you'll see below by the right multiply. The actual system is:

$$Ax + Bxv = c$$

At present, I am at a loss as to how to proceed, as the matrices aren't event the same dimensions. It is difficult for me to obtain the matrix $B$ easily, and is much easier and cheaper to compute matrix vector product like $Bc$, for the solution process.

EDIT: Before I said $A$ was $n\times n \times n$ and $B$ was $n\times n$, I got it switched.

  • $\begingroup$ Is it possible for you to calculate $A^{-1}$ or in another word is $A$ invertible or not? $\endgroup$ Dec 6, 2019 at 22:29
  • $\begingroup$ Oh sorry $A$ is not a square matrix :( $\endgroup$ Dec 6, 2019 at 22:32
  • $\begingroup$ Hey, sorry, A is square, but its large and non trivial to compute the inverse. I swapped the dimension of A and B by accident $\endgroup$
    – EMP
    Dec 6, 2019 at 22:52
  • $\begingroup$ What about an iterative method? Say find an initial guess by solving this: $Ax_{0} = c$ and then iteratively update by using this equation: $A x_{k+1} = c - B x_{k} v$ for $k > 0$? I have no idea it may work or not, so I will try to implement it. $\endgroup$ Dec 6, 2019 at 23:14
  • 2
    $\begingroup$ Have you tried to apply gmres to the whole system. You have a linear map $x \to Ax + Bxv$ and apply gmres on this. $\endgroup$
    – cfdlab
    Dec 7, 2019 at 5:46

3 Answers 3


Your second term is linear in $x$, so it can be rewritten as $Cx$ where $C$ is a suitable $n\times n$ matrix. You just need to figure out what its entries $C_{ij}$ are: for this you need a little index manipulation, but this should be an easy exercise in linear algebra.

Alternatively, you can use the (much better) notation $B(x\otimes v)$ for the second term (or $B(v\otimes x)$, depending on how your 'product' there is defined), and use Kronecker product properties to make the manipulations $(x\otimes v) = (I\otimes v)(x\otimes 1) = (I\otimes v)x$.

  • $\begingroup$ Thanks for the helpful answer. I unfortunately couldn't use it answer due to the issue I outlined in the question with only being able to compute $Bv$ rather than B outright. But I appreciate the clear response and I didn't make that issue clear enough. I accepted the answer as it did answer the question I asked (but I didn't make the implementation issue clear enough) $\endgroup$
    – EMP
    Dec 18, 2019 at 19:01

Based on Federico's answer, I think you can easily figure out what are $C_{ij}$ entries by explicitly writing the second term ($Bxv$) as:

$$(Bx)_{ij} = \sum_{\alpha=1}^{n} B_{ij\alpha} x_{\alpha}$$

$$(Bxv)_{i} = \sum_{\beta=1}^{n} (Bx)_{i\beta} v_{\beta}$$

$$(Bxv)_{i} = \sum_{\beta=1}^{n} \sum_{\alpha=1}^{n} B_{i\beta\alpha} x_{\alpha} v_{\beta} = \sum_{\alpha=1}^{n} \sum_{\beta=1}^{n} B_{i\beta\alpha}v_{\beta} x_{\alpha}$$


$$(Cx)_{i} = \sum_{\alpha=1}^{n} C_{i\alpha} x_{\alpha}$$


$$C_{i\alpha} = \sum_{\beta=1}^{n} B_{i\beta\alpha} v_{\beta}$$

So finally the main equation is written as:

$$Ax+Cx = c$$


$$(A+C)x = c$$

Now it might be possible to easily solve above linear equation by using direct or iterative methods.

  • $\begingroup$ Thanks for the clarification that you added to Federico Poloni's answer. I unfortunately couldn't use either answer due to the issue I outlined in the question with only being able to compute $Bv$ rather than B outright. But I appreciate the clear response and I didn't make that issue clear enough. I'm adding an answer explaining how I ended up tackling this that you can look to for the future. $\endgroup$
    – EMP
    Dec 18, 2019 at 19:01

So I got two great answers from Federico Poloni and Alone Programmer, but there was a problem in that all their solutions required having access to the full B matrix. Since in my case, the B matrix is the augmented hessian, or derivative of the augmented Jacobian used in BDF1-PTC nonlinear solvers, this is not a trivial matrix to obtain and I try to use the frechet derivatives to obtain the augmented Hessian product rather than the augmented hessian itself. As such I tried to find a way to do this with nested linear solves that would allow me to not have to compute the full Hessian and rely on the existing solver technology I have in my code.

One way to solve this could be by using a smoother like GS or Jacobi and solving the system: $$Dx = c - Bxv - Ax$$ Where D is the diagonal of the A matrix only. I chose not to employ this approach because GS does not have great solving properties. I instead chose to use my FGMRES with GS preconditioning like so:

  1. $x_0 = 0$
  2. Define $r = c - Bx_iv - Ax_i$
  3. Call FMGRES to solve $Ax_{i+1} = r$
  4. Define $r_2 = c - Bx_{i+1}v - Ax_{i+1}$
  5. If $r_2 < \epsilon $ stop, else i++ and go to 2

This iteration takes around 3 GMRES solves to fully converge the problem for $\epsilon = 1e-5$ and 6 for $\epsilon = 1e-15$. One thing to note is that early guesses don't need much precision. So an ideal implementation of this would have very low converging GMRES calls early on (or maybe use just some smoothing steps at first) before eventually using a more robust FGMRES with restart and many vectors for the final iterations.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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