Hot answers tagged

27

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 ...


21

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 ...


18

Use ArrayRules[]: m = SparseArray[{2, 2} -> i]; mc = SparseArray[ArrayRules[m] /. i -> -i, Dimensions[m]]; MatrixForm[mc] $\begin{pmatrix}0&0\\0&-\mathtt{i}\end{pmatrix}$


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 ...


17

I found a way to dramatically improve the performance of this algorithm by using the undocumented function SparseArray`KrylovLinearSolve. The key advantage of this function is that it seems to be a near-analog of MATLAB's pcg, and as such accepts as a first argument either: a square matrix, or a function generating a vector of length equal to the length ...


17

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, ...


16

J. M. has shown you a workaround using ArrayRules and as others mentioned, using Conjugate is more prudent. However, to answer your primary question — "Why doesn't ReplaceAll work on SparseArray?", it is because SparseArray is atomic. In other words, SparseArray objects are "indivisible" and the data contained in them can only be accessed in specific ways ...


16

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 :) *) ...


16

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@ ...


15

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] := ...


14

ArrayFlatten[{{sa11, sa12}, {sa21, sa22}}] seems to be what you need. It automatically merges everything into one big SparseArray[].


14

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 ...


13

array0 = {{72, 32, 64}, {18, 8, 16}, {63, 28, 56}}; array1 = SparseArray[Band[{# - 1 + Length@array0[[#]], #}, Automatic, {-1, 1}] -> array0[[#]] & /@ {1, 2, 3}, {5, 5}]; array1 // MatrixForm Update: Generalizing for arbitrary matrix input: rttF = Function[{mat}, With[{dims = Dimensions[mat]}, SparseArray[Band[{# - 1 + Last@dims, #}, ...


12

As Silvia notes in the comments, this behaviour is because of the predictive interface that is new in version 9. Mathematica tries to inspect the contents of the output to determine which contextual menu options to show in the suggestions bar depending on whether it is an array of integers/reals/mixed or if it is square/rectangular, etc. My guess is that ...


12

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 ...


12

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? ...


11

This question came up in Chat the other day. Here is the solution I proposed. banded[n_Integer?EvenQ] := With[ {main = RandomReal[99, n - 1], side = SparseArray[{}, n - 2, -0.5]}, SparseArray[{i_, i_} :> main[[i]], n - 1] + Sum[side ~DiagonalMatrix~ i, {i, {-1, 1}}] ] This uses several tricks and observations. Credit for the first ...


11

In addition to the answers given (ArrayFlatten), there may be cases when the SparseArray-s need to be accumulated one-by-one, and keeping them all in memory at the same time may be too expensive. I've written a generic function to gradually accumulate SparseArray-s. Here is the code: This is the low-level API for SparseArray construction / deconstruction, ...


11

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. ...


11

Here is a way using Band: SparseArray[{Band[{1, 2}, {4, 5}] -> {a, -a}, Band[{2, 1}, {5, 4}] -> {a, -a}}, {5, 5}] // MatrixForm $$\left( \begin{array}{ccccc} 0 & a & 0 & 0 & 0 \\ a & 0 & -a & 0 & 0 \\ 0 & -a & 0 & a & 0 \\ 0 & 0 & a & 0 & -a \\ 0 & 0 & 0 & ...


10

How about ArrayFlatten (docs): sa2 = ArrayFlatten[{{sa11, sa12}, {sa21, sa22}}] sa==sa2 (* True *)


10

This is a bug in Pick caused by SparseArray, has nothing to do with Graph. Minimal example (SparseArray object is the fullform version of your vLM): x = {1, 2, 3, 4, 5, 6}; Pick[x, SparseArray[Automatic, {6}, 0, {1, {{0, 3}, {{2}, {3}, {4}}}, {1, 1, 1}}], 0]; FullForm@x {1, System`Private`InternSequence[], System`Private`InternSequence[], ...


10

Re-applying SparseArray[] to a matrix or vector generated in this way usually restores the sparsity. p = SparseArray[{Band[{1, 1}] -> {1, 2, 4}, Band[{2, 1}] -> {5, -3}}, {3, 3}]; q = SparseArray[{{2, 2} -> -2, {3, 2} -> 3, {3, 1} -> -1}, {3, 3}]; r = p + q; rs = SparseArray[r]; Complement[r["NonzeroPositions"], rs["NonzeroPositions"]] ...


9

The determinant computation is a matter of memory use in terms of how much we want to store for subdeterminants of a Laplace expansion. Mathematica simply refuses to go that route after 11x11. YOu can do your own as below. myDet[mat_] /; Length[mat] <= 4 := Det[mat] myDet[mat_] := myDet[mat] = Sum[mat[[1, j]]*myDet[Drop[mat, {1}, {j}]], {j, ...


9

If you need a batch update, then the answer is in my comment you linked. If you need element-by-element, then there are two cases: Most of values you update are non-zero (or, generally, not equal to default element). In this case, I believe the answer of @Mr. Wizard is optimal, and you should expect update of a single element to be constant time. Most ...


9

Take a look at Unitize, UnitStep, Clip, Sign, etc. and the NonzeroPositions property of sparse arrays. Using these in combination with appropriate arithmetic operations you can quickly get positions. You could also take a look at FilterRules, and use this against the array rules to get positions for specified elements. With your updated post, e.g., ...


9

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.


9

You have shown that when the arrays are not sparse, using SparseArray is futile. Let's look at a case when it is sparse: MatSparse = SparseArray[{{1, 1} -> 1, {2000, 2} -> 2, {3, 3} -> 3, {1, 2000} -> 4}]; MatNormal = Normal[MatSparse]; multNormal = Timing[MatNormal.MatNormal]; multSparse = Timing[MatSparse.MatSparse]; Grid[{{"", "Normal", ...



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