58
Background
Details about multigrid solvers can be found in this pretty neat script by Volker John. That's basically the source from which I drew the information to implement the V-cycle solver below.
In a nutshell, a multigrid solver builds on two ingredients: A hierarchy of linear systems (with so-called prolongation operators mapping between them) and a ...
57
Introduction
This post is long overdue as I have been repeatedly asked to explain code of mine containing these things. As I see increased use of this construct by others perhaps it is past due also.
SparseArray objects can behave as functions accepting certain arguments to return internal data or efficiently return data in certain forms. These are known ...
33
Very good observation. Indeed, this issue is really frustrating. To single out the issue: It seems that Arnoldi's method is to blame:
Max@Abs@Im@Eigenvalues[M, -Nless]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Arnoldi"]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "Direct"]
Max@Abs@Im@Eigenvalues[M, -Nless, Method -> "FEAST"]
0.000610613
...
32
You can always modify the matrix so that the most negative eigenvalue is also the one with the largest absolute value, and hence corresponds to the first in the list returned by EigenVectors.
An upper bound for the largest absolute value of any eigenvalue is the Hilbert-Schmidt norm. So you can rescale your matrix by subtracting this norm times the unit ...
28
You are right, it can be done in a fraction of second. One can explicitly construct an array of indexes
blockArray[mat_] := SparseArray[
Tuples[Range@# - {1, 0, 0}].{Rest@#, {1, 0}, {0, 1}} &@Dimensions@mat ->
Flatten@mat]
Timings:
matrices = RandomReal[1, {48, 128, 128}];
s1 =
SparseArray@
ArrayFlatten@ReleaseHold@DiagonalMatrix[...
22
The way the code is written, you can neither exploit packed arrays nor any vectorization. There are two major reasons:
Using Rule prevents using packed arrays since Rule enforces symbolic computations. It is more efficient to feed SparseArray with pat -> vals, where pat is a packed array of dimensions {nnz, 2} of integers and vals is a packed array of ...
21
LinearSolve[] actually computes a permuted Cholesky decomposition; that is, it performs the decomposition $\mathbf P^\top\mathbf A\mathbf P=\mathbf G^\top\mathbf G$. To extract $\mathbf P$ and $\mathbf G$, we need to use some undocumented properties. Here's a demo:
mat = SparseArray[{Band[{2, 1}] -> -1., Band[{1, 1}] -> 2.,
Band[{1, ...
20
There is actually an undocumented System Option that tells Mathematica to do this automatically. The default behavior:
ind = {{3, 1}, {3, 3}, {1, 3}, {2, 1}, {3, 2}, {3, 1}, {3, 2}, {3, 3}, {1, 3}, {3, 1}};
val = {1, 1, 3, 0, 3, 4, 3, 1, 1, 1};
SparseArray[ind -> val] // Grid
$
\begin{matrix}
0 & 0 & 3 \\
0 & 0 & 0 \\
1 & 3 & ...
20
Using the properties of Block matrices:
$$\det\begin{pmatrix}\mathbf A&\mathbf B\\\mathbf C&\mathbf D\end{pmatrix}=\det(\mathbf A)\det\left(\mathbf D-\mathbf C\mathbf A^{-1}\mathbf B\right)$$
To visualize your matrix:
mat1 = mat;
{mat1[[;; 10, ;; 10]], mat1[[;; 10, 11 ;;]],
mat1[[11 ;;, ;; 10]], mat1[[11 ;;, 11 ;;]]} = Range@4; (* cool :) *)
Show[...
19
3D Example
The problem with direct solvers is that starting in 3 dimensions, their performance for dealing with matrices stemming from PDEs drops rapidly. This is why I wanted to show at least one 3-dimensional example. As there is no immediate analogon for Loop subdivision of tetrahedral meshes, I use hexahedral meshes instead.
Preparation
Here are some ...
18
Reposting my answer from here (its relevant part about SparseArray)
The anatomy of sparse arrays
We start with a generally useful API for construction and deconstruction of SparseArray objects:
ClearAll[spart, getIC, getJR, getSparseData, getDefaultElement, makeSparseArray];
HoldPattern[spart[SparseArray[s___], p_]] := {s}[[p]];
getIC[s_SparseArray] := ...
18
Yes, it is possible! There is a WorframGPULibrary (WGL), which I discover recently. It is undocumented, however there are beautiful examples in
$InstallationDirectory/SystemFiles/Links/CUDALink/CSource/
It is similar to LibraryLink, but allows CUDAMemory as an argument. I wrote the code below to call main CUSPARSE routines directly from Mathematica
$$
C :=...
18
In addition to Carl Woll's post:
Computing the pseudoinverse of a the graph Laplacian matrix (a.k.a. the KirchhoffMatrix) is very expensive and in general leads to a dense matrix that, if the graph is too large, cannot be stored in RAM. In the case that you have to compute only a comparatively small block of the resistance distance matrix, you can employ ...
17
You can use the first documented usage for SparseArray:
So what you want to do is collect all of the rules {i,j}->val before passing them to SparseArray. There is already a built-in method for efficiently collecting a set of rules with no duplicate keys, namely AssociateTo.
So with a slight modification of the OP's code:
data=<||>;
Do[
c = ...
16
Based on rcampion2012's answer to Efficient Implementation of Resistance Distance for graphs?, you could use:
resistanceGraph[g_] := With[{Γ = PseudoInverse[N @ KirchhoffMatrix[g]]},
Outer[Plus, Diagonal[Γ], Diagonal[Γ]] - Γ - Transpose[Γ]
]
Then, you can find the resistance using:
r = resistanceGraph[GridGraph[{10, 10}]];
r[[12, 68]]
1.60899
15
Pattern sparse arrays are used, for example, in NDSolve to specify the structure of a Jacobian. This can save significant time while integrating a differential equation. In other words pattern sparse arrays are useful if one knows something about the structure of a sparse array but not (yet) about the values. There are some examples in the FEM Programming ...
15
This is definitely faster. It took me some time to figure out the combinatorics and there might still be some potential for improvement within the compiled functions.
getColumnIndices = Compile[{{m, _Integer}, {n, _Integer}, {i, _Integer}},
Transpose[Partition[
Partition[
Join[
Table[k, {k, 1, i n}],
Flatten[Table[Table[k, {j, 1, ...
15
ClearAll[tminus1, tminus2, tminus3, tminus4, tminus5]
tminus1[n_] := ArrayPad[IdentityMatrix[n - 1, SparseArray], {{1, 0}, {0, 1}}]
tminus2[n_] := DiagonalMatrix[ConstantArray[1, n-1], -1]
tminus3[n_] := SparseArray[{i_, j_} /; j == (i - 1) -> 1, {n, n}]
tminus4[n_] := IdentityMatrix[n+1, SparseArray][[;;-2, 2;;]]
tminus5[n_] := Drop[IdentityMatrix[n+1, ...
14
The basic reason is that once you convert a tensor expression into a SparseArray, you've "given control" of all levels of that expression to SparseArray to manage on your behalf in an efficient way (the number of levels is the rank of the tensor, to mix jargon).
SparseArray will then try to maintain the illusion that those levels are still really there. ...
14
In version 10 Wolfram introduced some dynamic panels in output objects. These are intended to provide some useful brief information (I guess) about what is stored.
Because these panels are dynamic they trigger an unsafe dynamic content warning if the notebook is not in a trusted path. I prefer to switch this dynamic panelling off, which you can do by ...
14
List @@ s will do the trick:
s = 1000000 // RandomInteger[{1, 1000000}, {#, 2}] -> RandomReal[1, #]& // SparseArray;
s // Head
(* SparseArray *)
s // Dimensions
(* {1000000, 1000000} *)
l = List @@ s;
l // Head
(* List *)
l // Dimensions
(* {1000000, 1000000} *)
l // First // Head
(* SparseArray *)
Take[l, 4]
Why Does This Work?
Commentators ...
14
The fastest way to get the identity matrix as a sparse array is simply this:
IdentityMatrix[10000, SparseArray]; // AbsoluteTiming
Here are just some thoughts on why SparseArray is not the default for IdentityMatrix:
You don't always want SparseArray output when defining matrices, and it's not always possible to decide automatically whether a SparseArray ...
14
You cannot benefit from SparseArray if your matrix is not sparse! The data structure storing a SparseArray needs a certain amount of overhead which only pays off if your matrix is really, really sparse (density below a permille or so). Moreover, whenever you find yourself converting a large matrix to a sparse matrix you have actually done something wrong. ...
14
Be careful with assigning a "form" (e.g. MatrixForm) to a variable. Another observation is that a is a SparseArray uses rules to store the values and since ReplaceAll works with the FullForm of an expression care has to be taken (using Normal will make the array a regular matrix again).
This should work:
max = 1;
PotentialTilde[V0_] := SparseArray[{Band[{...
14
No, this is not built-in and I have been waiting quite long for such a feature.
PardisoLink
With the help of Szabolcs's package "LTemplate`" I had written such an interface for the Intel MKL Pardiso that is shipped with Mathematica a couple of years ago. I spent all day to bring it into a shape that has at least a small chance to run on a computer ...
13
Everything needs to be of the same precision including the background element.
Then this:
sp = SparseArray[{{2, 2} -> 1.}, {50, 50}, 0.];
ArrayFlatten[{{sp, 0.}, {0., -sp}}]
will return a SparseArray without creating a dense matrix first.
13
The complete list is (V11.2)
First or symbolic 0
Total or symbolic 1
List or symbolic 2
Anything else
For anything else:
SetSystemOptions[
"SparseArrayOptions" -> "TreatRepeatedEntries" -> Blah];
SparseArray[{{1, 2}, {1, 2}, {2, 2}} -> {2, 3, 4}] // Normal
{{0, Blah[2, 3]}, {0, 4}}
The reason these are undocumented is that it is not clear how ...
12
If the background value is 0 just run SparseArray on it:
sp = SparseArray[Band[{1, 1}] -> 1, {100, 100}];
ByteCount[sp]
(* 2152 *)
nrm = Normal[sp];
ByteCount[nrm]
(* 40168 *)
sp2 = SparseArray[nrm];
ByteCount[sp2]
(* 2152 *)
If it has another background value do SparseArray[nrm, Automatic, background], if you don't know the background value use the ...
12
Update: here Table is faster and more user-friendly then Array.
mat[n_] := LowerTriangularize@Table[2 (1 + Boole[j > 1]) (i - 1) Mod[i + j, 2],
{i, n}, {j, n}];
mat[10] // MatrixForm
It is fast and the result is packed array
mat[1000] // Developer`PackedArrayQ // AbsoluteTiming
(* {0.142522, True} *)
12
Not faster, but maybe easier to implement
makeKMat2[k1_, k2_] := Module[{A, m, n},
m = Length[k1];
n = Length[k2];
A = Map[SparseArray[Band[{1, 1}] -> #, {n, n}, 0.] &, k1, {2}];
With[{B = SparseArray@k2}, Do[A[[i, i]] += B, {i, 1, m}]];
ArrayFlatten[A]
]
Your implementation needs about 0.0048s; mine needs 0.0051s on my machine.
(b3m2a1 ...
Only top voted, non community-wiki answers of a minimum length are eligible