Matrix decomposition for a linear combination of matrices I want to solve the linear system,
$$(\mathbf{A}+\alpha \mathbf{B}) \mathbf{x} = \mathbf{b},$$
where $\mathbf{A,B}\in\mathbb{R}^{n\times n}$, $\mathbf{x,b}\in\mathbb{R}^{n\times 1}$ and $\alpha\in\mathbb{R}$.
Let say that I want to solve that linear system for different $\alpha$ values. Is there any $\mathbf{A}+\alpha \mathbf{B}$ decomposition that could be easily updated for different $\alpha$ values so I could solve fast each system? If not, is there any way to solve this set of systems fast?
PD: I do not need a full answer but just some hints in case it can be solved fast.
 A: Suppose that 
$$(A+\alpha B)x = b$$
$A,B \in \mathbb{R}^{n \times n} , x ,b  \in \mathbb{R}^{n} , \alpha \in \mathbb{R} $
I'm not sure how to preserve the original matrices if you just want the solution 
Let $C= A+\alpha B$ 
now $C =U \Sigma V^{T} $
$$ Cx=b \implies U\Sigma V^{T}x =b$$
$$  C^{+} = V \Sigma^{-1} U^{T}  $$
$$ x = V \Sigma^{-1}U^{T} b $$
You should understand that even slight changes in $\alpha$ could change the condition number $(A+\alpha B)$ for some matrices which would destroy your precision. The computational time for solving this is very little to update it.
A: This is a very partial answer.
First, I am working over an arbitrary algebraically closed field $\mathbb{K}$.  Let $\mathbf{A}$ and $\mathbf{B}$ be $n$-by-$n$ matrices over $\mathbb{K}$.  The only situation in which you can invert $$\mathbf{X}(\alpha):=\mathbf{A}+\alpha\,\mathbf{B}$$
for any $\alpha \in\mathbb{K}$ is when $\mathbf{A}$ is invertible and $\mathbf{A}^{-1}\,\mathbf{B}$ is nilpotent.
Let $\mathbf{b}\in\mathbb{K}^n$ be given.  Thus, the only $\mathbf{x}(\alpha)\in\mathbb{K}^n$ satisfying $$\mathbf{X}(\alpha)\,\mathbf{x}(\alpha)=\mathbf{b}$$ must satisfy
$$\left(\mathbf{I}+\alpha \,\mathbf{A}^{-1}\,\mathbf{B}\right)\,\mathbf{x}(\alpha)=\mathbf{A}^{-1}\,\mathbf{b}\,,$$
where $\mathbf{I}$ is the $n$-by-$n$ identity matrix.
Because $\mathbf{A}^{-1}\,\mathbf{B}$ is a nilpotent $n$-by-$n$ matrix, we conclude that $\left(\mathbf{A}^{-1}\,\mathbf{B}\right)^n$ is the zero matrix.  Ergo,
$$\mathbf{x}(\alpha)=\big(\mathbf{I}+\alpha\,\mathbf{A}^{-1}\,\mathbf{B}\big)\,\mathbf{A}^{-1}\,\mathbf{b}=\sum_{k=0}^{n-1}\,(-1)^k\,\alpha^k\,\left(\mathbf{A}^{-1}\,\mathbf{B}\right)^{k}\,\mathbf{A}^{-1}\,\mathbf{b}\,,$$
where we use the convention $t^0:=1$ and $\mathbf{M}^0:=\mathbf{I}$ for any $t\in\mathbb{K}$ and any $n$-by-$n$ matrix $\mathbf{M}$ over $\mathbb{K}$.  Thus, the only lengthy computation is the computation of $\mathbf{A}^{-1}$ (and $\mathbf{A}^{-1}\,\mathbf{B}$).

If you work over $\mathbb{R}$ or $\mathbb{C}$, then you can fix an operator norm $\|\_\|$.  The easiest choice is perhaps the Frobenius norm.  When $\mathbf{A}$ is invertible and $$|\alpha|\,\left\|\mathbf{A}^{-1}\,\mathbf{B}\right\|<1\,,$$ then $\mathbf{X}(\alpha)$ is invertible and you can approximate $\mathbf{x}(\alpha)$ using the infinite sum below:
$$\mathbf{x}(\alpha)=\sum_{k=0}^\infty\,(-1)^k\,\alpha^k\,\left(\mathbf{A}^{-1}\,\mathbf{B}\right)^k\,\mathbf{A}^{-1}\,\mathbf{b}\,.\tag{*}$$
The convergence will be painfully slow when $\delta:=|\alpha |\,\left\|\mathbf{A}^{-1}\,\mathbf{B}\right\|$ is very close to $1$.  Nevertheless, if you want a certain amount of accuracy (say, within $\epsilon$), you can pre-compute $\left(\mathbf{A}^{-1}\,\mathbf{B}\right)^k$ up to $$k\approx \dfrac{\ln(1-\delta)+\ln(\epsilon)}{\ln(\delta)}$$ and be done with the job for all $\alpha$ with $|\alpha|\leq \dfrac{\delta}{\left\|\mathbf{A}^{-1}\,\mathbf{B}\right\|}$.

The field $\mathbb{K}$ is either $\mathbb{R}$ or $\mathbb{C}$.  For some special $\alpha$, we can use a similar power-series approach as in (*).
If $\mathbf{A}$ is singular but $\mathbf{B}$ is invertible, then you can just look at
$$\mathbf{Y}(\beta):=\mathbf{B}+\beta\,\mathbf{A}\,.$$
For $\alpha\neq 0$, we see that
$$\alpha\,\mathbf{Y}\left(\frac{1}\alpha\right)=\mathbf{X}(\alpha)$$
and we can do a similar trick to (*) with $\mathbf{Y}(\beta)$ instead.
In the case where both $\mathbf{A}$ and $\mathbf{B}$ are singular, you simply have to find an element $\lambda\in\mathbb{K}$ such that $\mathbf{C}:=\mathbf{A}+\lambda\,\mathbf{B}$ is nonsingular.  If $\lambda$ does not exist, then you will be solving degenerate linear system for all $\lambda$, which is obviously no fun.  Nonetheless, if you are lucky, this $\lambda$ exists, and we can use
$$\mathbf{Z}(\gamma):=\mathbf{C}+\gamma\,\mathbf{B}$$
instead, noting that $\mathbf{X}(\alpha)=\mathbf{Z}(\alpha-\lambda)$.  Then, you can work with $\mathbf{Z}(\gamma)$ just like before.  
