11

Absolutely! First off, some linear algebra systems are smart enough to only store half of the matrix, this could save you a bunch of memory. But even if this isn't the case, various algorithms in numerical linear algebra will exploit the symmetry. For example, given a symmetric matrix, any eigensolver will immediately know that all eigenvalues are real-...


9

What one usually wants in this situation is to preserve a discrete analog of time symmetry: namely, if the time discretization is applied to solve first forward and then backward in time, the initial condition is recovered. This is true if the method is invariant under the following substitutions: $$ \Delta t \to -\Delta t$$ $$ a^{n+j} \to a^{n-j}$$ (here ...


8

Let $P$ be the anti-diagonal permutation matrix, $$P = \begin{bmatrix} & & & 1 \\ & & 1 \\ & 1 \\ 1 \end{bmatrix}$$ so that $PAP$ is the version of $A$ with rows and columns reversed. The first $P$ swaps the rows, and the last $P$ swaps the columns. We have the Cholesky decomposition $$PAP=LL^T$$ which implies $$A = (PLP)(PLP)^T,$$ ...


6

Yes, for an SPD matrix $\mathbf A$ there are a variety of Cholesky-like decompositions, you can derive the $\mathbf A = \mathbf L^T \mathbf L$ variant by first writing down an educated/structured guess.. $\begin{bmatrix} \mathbf A_{11} & \mathbf a_{21}^T \\ \mathbf a_{21} & \alpha_{22} \end{bmatrix} = \begin{bmatrix} \mathbf L_{11}^T & \...


6

The CG method works by producing in every iteration a vector $x^k\in\mathbb{R}^n$ that solves the minimization problem (assuming $x^0=0$ for simplicity) $$\min_{x\in \mathcal{K}_k} \frac12(x-x^*)^TA(x-x^*),\tag{1}\label{cg1}$$ where $x^*=A^{-1}b$ and $\mathcal{K}_k$ is a Krylov space of (apart from exceptional circumstances) dimension $k$. (The CG algorithm ...


6

We have the matrix $A$ that can be expressed as $A = JGJ^T$. The first thing is to calculate the QR decomposition of matrix $J$. Because of the low rank of the matrix it can be done very fast with, for instance, modified Gram Schmidt algorithm. Now we can write $A$ as $A = QR G R^TQ^T$, where $Q$ is an orthonormal matrix ($Q^T Q = I$). We define $F$ as ...


5

As the comments note, it shouldn't be an issue with your eigensolver if the matrix is symmetric to machine precision. That being said, however, if you are really worried about it, just copy the transpose of the upper triangle to lower triangle after you generate the mass matrix, and then it will be exact.


4

Since a SPD matrix is invertible, we can make the Cholesky decomposition $A^{-1} = PP^T$. Since $A$ is non-singular, so is $P$, and the inverse of a triangular matrix is triangular, so writing $L = P^{-1}$ we have $A^{-1}$ = $L^{-1}L^{-T}$. Inverting both sides gives $A = L^TL$.


4

In short, orthogonalization of the Krylov vectors occurs with respect to the operator, but not with respect to the preconditioner. Alright, so say we want to solve $Ax=b$ with preconditioner $B$. the preconditioned-CG iteration is basically: \begin{align*} \hat{v}_1=\tilde{v}_1 =& Bb\\ v_1 =& \tilde{v}_1 / c_1\\ \\ \hat{v}_i =& BAv_{i-1}\\ \...


2

This is a typical case where it does not make sense to actually store the matrix. Rather, consider the matrix to be an operator where if you need to multiply by it, i.e., form $y=Ax$, you do $y=b (b\cdot x)$ instead if you need to access an element $A_{ij}$, you instead compute it as $b_ib_j$. In other words, knowing what $A$ is, you should just not store ...


2

The differential operator is "grad after div" and so (in the continuous version) is symmetric and negative semi-definite. To get a symmetric discretization your discretization of the gradient has to be the (negative) adjoint of your discrete divergence. With finite differences this says, for example, that forward differences for the divergence imply backward ...


2

It would seem to me that you can subdivide your domain into a relatively small set of tetrahedra. Then, it is trivial to do the integration because there are many good (and pre-tabulated) quadrature rules on tetrahedra that will yield reasonably high accuracy -- this is what one does in the finite element method all the time, so there is a lot of information ...


1

You may want a method that works right down to the coordinate singularity at $r=0$. I will do the spherical case, but the cylindrical case is similar. We want to avoid ever dividing by $r$. To think of a Finite Volume method, consider the control volumes to be concentric shells. The inner and outer radii of shell i are $(r_i-\Delta r/2, r_i+\Delta r/2)$, ...


1

I will assume you mean complex symmetric s.t. $A = A^T$ and not Hermitian $A = A^H$, which is different. Let $B \in \mathbb{R},[2n \times 2n]$ be your new tri-diagonal matrix. A simple option is each imaginary component is mapped $n$ places down the appropriate diagonal: Main diag.: $b_{i,i} = Re(a_{i,i}), \;b_{i+n,i+n} = Im(a_{i,i}) \qquad \text{for} \...


1

Read http://www.sciencedirect.com/science/article/pii/S1877050915010492 "Nonsymmetric Preconditioning for Conjugate Gradient and Steepest Descent Methods"


1

It is possible to use any group $G$ for isomorphism pruning with Margot's isomorphism pruning (ISOP 1.1) solver. One way of doing this is by creating a different integer linear program (ILP) whose symmetry group is $G$ and submitting this ILP to Margot's solver for finding the symmetry group. Then you can replace this ILP with the original ILP and submit ...


1

I think you should explore the geometry of the object described by your permutation group. Consider for example if you had three variables and if your permutations allow to permute $x_2$ and $x_3$. Then your optimization problem is posed on the hyperplane $(x_1,x_2,x_2)$. Similarly, if you allowed the full permutation group, you'd be optimizing on the line $...


1

At this point quite a few papers have been written about exploiting symmetry in integer programming and symmetry breaking techniques have been implemented by lots of people and are available in the widely used CPLEX and Gurobi solvers. So yes, these techniques can be useful in practice. Your question is really quite vague. Could you be more specific ...


1

Can you use timestepping instead of discretizing in time? You have a second order time derivative, but you can rewrite that as a first order system, i.e. $\frac{\partial^2 u}{\partial t^2} = f$ can be rewritten by adding another equation $\frac{\partial w}{\partial t} = f$ $\frac{\partial u}{\partial t} - w = 0,$ and then use time-stepping methods for ...


1

Many plane-wave codes use a very simple k-point weighting scheme: Generate a uniform mesh in k-space, $K$ Assign each k-point in the irredicuble wedge a weight $w(k) = \#\{q \in K : \text{k and q are symmetric}\}$ Write integrals as $\int_{BZ} F(k) dk = \sum_{k \in K} w(k) F(k)$ This is quick and dirty, but it usually isn't the leading error (for PWDFT, at ...


1

An LDLT decomposition will give you the information you need. If the diagonal $D$ matrix from the LDLT decomposition of your matrix is positive definite, your matrix is also positive definite. Here is a link to some documentation discussing the LDLT decomposition in Eigen: http://eigen.tuxfamily.org/dox-2.0/TutorialAdvancedLinearAlgebra.html


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