We’re rewarding the question askers & reputations are being recalculated! Read more.
3 corrected term
source | link

Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace equationoperator, $\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 example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

julia> using SpecialMatrices
julia> Strang(5)
5×5 SpecialMatrices.Strang{Float64}:
 2.0  -1.0   0.0   0.0   0.0
-1.0   2.0  -1.0   0.0   0.0
 0.0  -1.0   2.0  -1.0   0.0
 0.0   0.0  -1.0   2.0  -1.0
 0.0   0.0   0.0  -1.0   2.0

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "pseudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

However, let's say we want to invert the matrix.

julia> inv(collect(Strang(5)))
5×5 Array{Float64,2}:
 0.833333  0.666667  0.5  0.333333  0.166667
 0.666667  1.33333   1.0  0.666667  0.333333
 0.5       1.0       1.5  1.0       0.5
 0.333333  0.666667  1.0  1.33333   0.666667
 0.166667  0.333333  0.5  0.666667  0.833333

Notice that the inverse of this sparse matrix is dense. Thus, if we do not know the analytical solution of the inverse (which is true for most sparse matrices arising from PDE discretizations), the amount of memory that will be required for the inverse is $\mathcal{O}(n^2)$. This means that the inverse will take up a massive amount of extra memory, and your memory limit will thus be determined by the size of the dense inverse matrix.

Instead, direct solver methods of \ and iterative solvers like those in IterativeSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

As others have mentioned, condition number and numerical error is another reason, but the fact that the inverse of a sparse matrix is dense gives a very clear "this is a bad idea".

Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace equation, $\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 example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

julia> using SpecialMatrices
julia> Strang(5)
5×5 SpecialMatrices.Strang{Float64}:
 2.0  -1.0   0.0   0.0   0.0
-1.0   2.0  -1.0   0.0   0.0
 0.0  -1.0   2.0  -1.0   0.0
 0.0   0.0  -1.0   2.0  -1.0
 0.0   0.0   0.0  -1.0   2.0

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "pseudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

However, let's say we want to invert the matrix.

julia> inv(collect(Strang(5)))
5×5 Array{Float64,2}:
 0.833333  0.666667  0.5  0.333333  0.166667
 0.666667  1.33333   1.0  0.666667  0.333333
 0.5       1.0       1.5  1.0       0.5
 0.333333  0.666667  1.0  1.33333   0.666667
 0.166667  0.333333  0.5  0.666667  0.833333

Notice that the inverse of this sparse matrix is dense. Thus, if we do not know the analytical solution of the inverse (which is true for most sparse matrices arising from PDE discretizations), the amount of memory that will be required for the inverse is $\mathcal{O}(n^2)$. This means that the inverse will take up a massive amount of extra memory, and your memory limit will thus be determined by the size of the dense inverse matrix.

Instead, direct solver methods of \ and iterative solvers like those in IterativeSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

As others have mentioned, condition number and numerical error is another reason, but the fact that the inverse of a sparse matrix is dense gives a very clear "this is a bad idea".

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 example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

julia> using SpecialMatrices
julia> Strang(5)
5×5 SpecialMatrices.Strang{Float64}:
 2.0  -1.0   0.0   0.0   0.0
-1.0   2.0  -1.0   0.0   0.0
 0.0  -1.0   2.0  -1.0   0.0
 0.0   0.0  -1.0   2.0  -1.0
 0.0   0.0   0.0  -1.0   2.0

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "pseudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

However, let's say we want to invert the matrix.

julia> inv(collect(Strang(5)))
5×5 Array{Float64,2}:
 0.833333  0.666667  0.5  0.333333  0.166667
 0.666667  1.33333   1.0  0.666667  0.333333
 0.5       1.0       1.5  1.0       0.5
 0.333333  0.666667  1.0  1.33333   0.666667
 0.166667  0.333333  0.5  0.666667  0.833333

Notice that the inverse of this sparse matrix is dense. Thus, if we do not know the analytical solution of the inverse (which is true for most sparse matrices arising from PDE discretizations), the amount of memory that will be required for the inverse is $\mathcal{O}(n^2)$. This means that the inverse will take up a massive amount of extra memory, and your memory limit will thus be determined by the size of the dense inverse matrix.

Instead, direct solver methods of \ and iterative solvers like those in IterativeSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

As others have mentioned, condition number and numerical error is another reason, but the fact that the inverse of a sparse matrix is dense gives a very clear "this is a bad idea".

2 Improve writing
source | link

Here'sHere is a quick example which is very practical related to memory usage in PDEs. When one discretizes a LaPlacian $ \Delta u $ saythe Laplace equation, $\Delta u$, for example, in the Heat Equation

$$ u_t = \Delta u + f(t,u) $$$$ u_t = \Delta u + f(t,u) .$$

in order toTo solve it numerically, one ends up with sparse matrices $A$, and a MOLmethod of lines discretization then solve

The canonical 1D example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jlSpecialMatrices.jl:

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "psudo"pseudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

Instead, direct solver methods of \ and iterative solvers like those in IteartiveSolvers.jlIterativeSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

Here's a quick example which is very practical related to memory usage in PDEs. When one discretizes a LaPlacian $ \Delta u $ say in the Heat Equation

$$ u_t = \Delta u + f(t,u) $$

in order to solve it numerically, one ends up with sparse matrices $A$, and a MOL discretization then solve

The canonical 1D example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "psudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

Instead, direct solver methods of \ and iterative solvers like those in IteartiveSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace equation, $\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

The canonical 1D example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "pseudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

Instead, direct solver methods of \ and iterative solvers like those in IterativeSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

1
source | link

Here's a quick example which is very practical related to memory usage in PDEs. When one discretizes a LaPlacian $ \Delta u $ say in the Heat Equation

$$ u_t = \Delta u + f(t,u) $$

in order to solve it numerically, one ends up with sparse matrices $A$, and a MOL discretization then solve

$$ u_t = Au + f(t,u) $$

The canonical 1D example is the Strang matrix. Implicit methods will need to invert $A$, or some form of $I-\gamma A$. For a discretization with 5 points in space, let's see what happens to this operator. We can generate it easily in Julia using SpecialMatrices.jl:

julia> using SpecialMatrices
julia> Strang(5)
5×5 SpecialMatrices.Strang{Float64}:
 2.0  -1.0   0.0   0.0   0.0
-1.0   2.0  -1.0   0.0   0.0
 0.0  -1.0   2.0  -1.0   0.0
 0.0   0.0  -1.0   2.0  -1.0
 0.0   0.0   0.0  -1.0   2.0

This special matrix is tridiagonal (many other discretizations are banded, and thus still sparse). This means that you can store it by only storing 3 arrays, or in this case, it can be "lazy" (i.e., no array is needed and you can have a "psudo-array" type which generates the values on demand). This means that even for very large spatial discretizations with $n$ points, sparse matrix formats can store this in $\mathcal{O}(3n)$ memory (and lazy formats can do this in $\mathcal{O}(1)$!).

However, let's say we want to invert the matrix.

julia> inv(collect(Strang(5)))
5×5 Array{Float64,2}:
 0.833333  0.666667  0.5  0.333333  0.166667
 0.666667  1.33333   1.0  0.666667  0.333333
 0.5       1.0       1.5  1.0       0.5
 0.333333  0.666667  1.0  1.33333   0.666667
 0.166667  0.333333  0.5  0.666667  0.833333

Notice that the inverse of this sparse matrix is dense. Thus, if we do not know the analytical solution of the inverse (which is true for most sparse matrices arising from PDE discretizations), the amount of memory that will be required for the inverse is $\mathcal{O}(n^2)$. This means that the inverse will take up a massive amount of extra memory, and your memory limit will thus be determined by the size of the dense inverse matrix.

Instead, direct solver methods of \ and iterative solvers like those in IteartiveSolvers.jl will solve $Ax=b$ without ever computing $A^{-1}$, and thus only have the memory requirement of having to store $A$. This can greatly expand the size of the PDEs you can solve.

As others have mentioned, condition number and numerical error is another reason, but the fact that the inverse of a sparse matrix is dense gives a very clear "this is a bad idea".