60

Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops. $v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a ...


19

Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table: operation memory reads/writes matrix + matrix 3n² matrix * vector 2n²+n (if vector is not cached) matrix * vector n²+2n (if vector is only read once) vector + vector ...


4

I think you already know the answer- find the inverse of each of the small diagonal blocks. The details of how you might break the matrix up into the block in Python are something out of scope for this stackexchange group. As others have suggested you most likely don't really need the inverse of this matrix and could probably do everything you need with a ...


4

You might want to use the fact that: $$ ||A||_2=\sigma_\max(A) $$ where $\sigma_\max$ is the largest singular value. If you are interested in details, this Math SO question should be interesting. Thus, $$ ||A^{-1}||_2=\frac{1}{\sigma_\min(A)} $$ where $\sigma_\min$ is the smallest singular value. You certainly want to avoid the actual calculation of the ...


3

For me, a classic problem in linear algebra that got me climbing down a rabbit hole of complex analysis, conformal mappings, and polynomials, was in trying to prove converge rate bounds and iteration bounds for Krylov subspace methods. At a high level, it asks one of the following questions: What is the fastest way of computing $A^{-1}b$, under the ...


3

Short answer: When he says that when $x=0$, any $u$ will do, he means any $u$ that satisfies $\|u\|_2 = \sqrt{2}$. Otherwise, it wouldn't be a Householder reflector. So yes, you have found a transformation that "introduces" zeros in $x$, but you haven't proven that that transformation is a Householder reflector (spoiler: it isn't ;) ).


2

By definition, the Householder transformation is a reflection about a plane or hyperplane. The plane is described by its unit normal vector u and the transformation is then H = I - 2uu'. If u is a unit vector, 2uu' cannot be zero and thus H cannot be the identity matrix. You can find plenty of non-Householder transformations H that, for certain arguments, ...


2

As I mentioned in my comment, due to that you are searching for a method based on Python and ideally available in NumPy or SciPy, my suggestion is to use scipy.linalg.expm. It's really easy to use basically you have the $X$ matrix and you just pass it to the scipy.linalg.expm class and it would give you the exponential of $X$ (i.e. $e^{X}$). Due to that you ...


2

If changes in pivoting are an issue, then yes, to my knowledge there is no obvious $O(n^3)$ solution. However, you could consider switching to the QR factorization. That factorization costs twice as much as LU for the initial computation, but it can be updated in a stable way in $O(n^2)$. See Golub and Van Loan Matrix Computations for algorithms, or https://...


2

The Newton-Raphson method can be used to solve non-linear systems of equations. The first step is to write your system as a root finding problem: $$ f(T_n) = \left( \frac{C}{\Delta t} + K \right) T_n - Q_n - \frac{C}{\Delta t} T_{n-1} = 0 $$ Taylor expand this equation about an initial guess $T^0_n$, keeping only the linear term: $$ f(T_n) \approx f(T^0_n)...


2

First of all, I find the name "Smoothed Aggregation" a bit misleading, because the method - as I understood it - consists of both smoothing a tentative prolongation operator and implicitly considering the (near-)null-space for constructing the tentative prolongation operator. In standard algebraic multigrid with linear elasticity, corrections from coarser ...


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