Hot answers tagged

19

(This is getting too long for comments...) I'll assume you actually need to compute an inverse in your algorithm.1 First, it is important to note that these alternative algorithms are not actually claimed to be faster, just that they have better asymptotic complexity (meaning the required number of elementary operations grows more slowly). In fact, in ...


12

Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace operator, $\Delta u$, for example, in the Heat Equation $$ u_t = \Delta u + f(t,u) .$$ To solve it numerically, one ends up with sparse matrices $A$, and a method of lines discretization then solve $$ u_t = Au + f(t,u) $$ The canonical 1D ...


12

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


11

Since $$ A = B(I-B)^{-1} = (I-B)^{-1}(I-B)B(I-B)^{-1} = (I-B)^{-1}B(I-B)(I-B)^{-1} =(I-B)^{-1}B $$ So you want to solve $$ (I-B)A=B $$ You seem to need only the first three columns of $A$. Solve the matrix problems $$ (I-B)a_i = b_i, \qquad i=0,1,2 $$ where $b_0,b_1,b_2$ are first three columns of $B$. Then $a_0,a_1,a_2$ are the first three columns of $A$.


11

Normally there are some principal reasons to prefer solve a linear system respect to use the inverse. Briefly: problem with the conditional number (@GoHokies comment) problem in the sparse case (@ChrisRackauckas answer) efficiency (@Kirill comment) Anyway, as @ChristianClason remarked in comments, can be some cases where the use of the inverse is a good ...


8

Your two ideas make it much too complicated. If $X$ is the inverse of $A$, $$ AX=I, $$ and $x_i$ is the $i$-th column of $X$ and $e_i$ is the $i$-th column of the identity matrix $I$ ($e_i$ is a vector of all zeros except with $1$ in the $i$-th location), then the columns $x_i$ of the inverse are defined by $$ Ax_i = e_i. $$ All you need to do is solve $n$ ...


7

Think of your matrix as block-diagonal with blocks $$ C_\ell = \begin{pmatrix} 0 & & \\ & D_\ell & \\ & & 0\end{pmatrix}. $$ Then it is clear that $$ C_\ell^{-1} = \begin{pmatrix} 0^{-1} & & \\ & D_\ell^{-1} & \\ & & 0^{-1}\end{pmatrix} $$ where of course $0^{-1}$ is not well defined. But that doesn't ...


6

As suggested in the comments, the eigendecomposition $\mathbf J = \mathbf Q \mathbf D \mathbf Q^{T}$ can be used to generate the matrix $\mathbf J^{-1/2}$, just take the eigenvectors from $\mathbf J$ (denoted $\mathbf Q$) and the inverse square root of the eigenvalues (denoted $\mathbf D^{-1/2}$) .. this is a "scalar"/"entrywise" operation because $\mathbf ...


6

Pseudoinverses typically will be computed via some truncation procedure to determine the rank, so they are not close to the original inverse. Example: $$A = Q \begin{bmatrix} 1\\ & 10^{-12} \\ & & 10^{-18} \end{bmatrix}Q^{*} $$ has exact (pseudo)inverse $$A^+ = A^{-1} = Q \begin{bmatrix} 1\\ & 10^{12} \\ & & 10^{18} \end{bmatrix}Q^{*},...


5

I will try to give my thought on the first question regarding fast $3\times 3$ inverse. Consider $$ A=\left[ \begin{array}{ccc} a & d & g\\ b & e & h\\ c & f & i \end{array}\right] $$ Since the matrices are small and very general (do not feature any known structure, zeroes, relative scales of the elements), I think it would be ...


5

Unfortunately, although it is easy to update the Cholesky factorization after a low-rank update, updating the Cholesky factorization after a full rank update like this isn't any easier than recomputing the Cholesky factorization. An alternative to consider would be to compute the eigenvalue decomposition of your matrix, $A=LL^{T}=UDU^{T}$. Adding $\mu^{...


5

I think this is a really cool question, which might have a really cool answer, if someone was willing to think about it hard enough to publish on. But when we're talking smooth functions which need to be evaluated quickly, I reach for B-splines. Create the pairs $\{x_j, f(x_j)\}_{j=1}^{n}$, and use (say) a cubic B-spline (here's an implementation for the ...


4

You should probably note that, buried deep inside the numpy source code (see https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src) the inv routine attempts to call the dgetrf function from your system LAPACK package, which then performs an LU decomposition of your original matrix. This is morally equivalent to Gaussian elimination, but ...


4

The Cholesky decomposition [the function dpotrf() in LAPACK] factors $\mathbf A = \mathbf L \mathbf L^{\mathrm T}$, or alternatively $\mathbf A^{-1} = \left(\mathbf L \mathbf L^\mathrm T \right)^{-1} = \mathbf L^{-\mathrm T}\mathbf L^{-1}$. If you insert the latter representation your other expressions you'll see how you can compute them efficiently: $\...


3

So, a cofactor matrix is a transpose of an adjugate matrix. I know of the following paper: G. W. Stewart, "On the adjugate matrix," Lin. Alg. Appl., vol. 283, no. 1–3, pp. 151–164, Nov. 1998. (available with full text) There, the author works on an algorithm of computing a adjugate matrix $\text{adj}(A)$ when $A$ is nearly singular or singular. For such ...


3

Depending on the condition number of the matrix $A$, computing Padé approximations of the matrix exponential $e^A$ is horribly ill-conditioned, especially for such high orders. If a problem is ill-conditioned, no amount of mathematical trickery will allow a stable computation; however, often you can change the problem to one that is better conditioned (and ...


3

You might be better served by either the LDL' decomposition or the Cholesky decomposition (in the event that C's are positive definite in addition to symmetric, they probably are). Though all of the algorithms have cubic complexity, the constants are much better for these "linear" solvers than the SVD (which is more of a "spectral" solver). They are also ...


3

The pseudo inverse $A^+$ fulfills: $A^+b$ is the minimimum norm solution of the least squares problem $\min_x \|Ax-b\|_2^2$. Hence, to calculate $A^+e_k$ by solving the respective optimization problem (which has a different structure, depending on the format of $A$, as Brian Borchers commented. You could also use an iterative method to solve $A^TA x = A^...


3

There are a number of symbolic algebra packages that will compute the exact matrix inverse; e.g. Sage, Mathematica, Maple, or Maxima. This becomes computationally expensive for large matrices, but there is a chance that a sparse $100 \times 100$ matrix would be within reach. I wouldn't rule it out without trying, since it's so easy to try.


3

You don't need the inverse even with the goal of finding $K^{-1} h h^{T} K^{-1} - K^{-1}$. If you are interested to have this expression, I would explain how you can convert it to a matrix equation and then solve it more efficiently: Let's define the $X$ as: $$X = K^{-1} h h^{T} K^{-1} - K^{-1}$$ Your objective is to calculate $X$ in this equation by ...


2

This may be more of a comment than an answer, but I can't comment. Yes on arbitrary precision. No on symbolic inverse, because at some point, you have to numerically evaluate it, and it is likely to be extremely numerically unstable. Even for 6 by 6 matrix, the symbolic inverse is starting to get pretty ungainly. Which leaves me the main item. Why do you "...


2

Perhaps it would help to decompose $v$ in terms of the (normalized) eigenvectors $w_i$ of the matrix $A$. Let the corresponding eigenvalues be $\lambda_i$, so that $A \cdot w_i = \lambda_i w_i$. If $c_i = <w_i, v>$ then $(I + A + A^2 + ... + A^n) \cdot v = \sum_i c_i (1+\lambda_i + \lambda_i^2 + ... + \lambda_i^n) w_i$ For the eigenvalues that have $|...


2

Unfortunately, the listed set of matrix properties does not give any hope on a closed-form of the matrix inverse or its reasonable approximation (for this specific case). It is too general, and finding matrix inverse (or matrix factorization) is a non-trivial task. Otherwise, such fundamental problem as preconditioning would have been successfully solved. ...


2

You can rephrase your question as an optimization problem and solve numerically with gradient descent. You would probably converge quite fast, because your initial guess will be the solution you already have for $\mu^2=0$. More formally, we like to find $x$ that minimizes: $f(x) = ||(LL^T+μ^2I)x -y||^2$, where $L, μ^2, y$ are given. To solve it, start ...


2

One easy way to derive finite-difference approximations to derivatives is as follows. To find the coefficients $c_k$ corresponding to order-$s$ derivative on points $uh,(u+1)h,\ldots,vh$ ($u=-10,v=10,s=1$ in your case), write the condition defining, using Taylor theorem, as $$ \sum_{j=u}^{v}c_k e^{j h \partial} = \partial^s, $$ which should hold ...


2

Cramer's rule using a more stable determinant: https://hal.archives-ouvertes.fr/hal-01500199/document Implementation of the algorithm is consistently 2.5 times slower than LU factorization of matlab regardless of size of matrix.


1

I known that this is not a definitive answer, but it can give an idea how to move with CUDA. At this stage is difficult to give advice because for this is necessary a detailed know about the actually code As already write in comments, in general, is better not invert the matrix, but solve the linear system (for detail about why see this question). There is ...


1

Three things immediately come to mind: R might not take advantage of sparsity when using the solve command to compute the inverse of a matrix. Usually, the inverse of a sparse matrix is dense, so the built-in procedure might have been written to just convert it to a dense matrix, figuring that there's rarely anything to be gained from sparsity. To check ...


1

The problem is wildly underdetermined. Of the infinitely many equations for $J$ of the equation $J^T e=f$ is the (rank-one) matrix $$ J^T = \frac{f e^T}{\|e\|^2}. $$ For this, you have that $$ JJ^T = \frac{e f^Tf e^T}{\|e\|^4} = \|f\|^2 \frac{e e^T}{\|e\|^4}. $$ Interestingly, this matrix satisfies $$ (JJ^T)(JJ^T) = \|f\|^4 \frac{e e^T}{\|e\|^4} \frac{...


1

I would not generally expect a "20th order" derivative estimate to typically be very stable/reliable/useful (e.g. due to well known artifacts of high-order polynomial interpolation). That said, a general procedure for deriving finite-difference stencils is to solve an appropriate polynomial interpolation problem. I will use Matlab-style notation, as you ...


Only top voted, non community-wiki answers of a minimum length are eligible