Hot answers tagged

18

Here's an edited version of my answer to a related question (elsewhere). Since your central question was about speed (or time complexity), you might wish to know an important result from elementary theory of algorithms and computational complexity, which is that the time and space complexity of matrix multiplication depends upon the order of such ...


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

This happens because of unpacking when the numbers exceed $MaxMachineNumber: fast = Dot @@@ Partition[tab, Divisors[3960][[42]]]; Developer`PackedArrayQ /@ fast (* {True, True, True, True, True, True, True, True} *) Max[fast] <= $MaxMachineNumber (* True *) slow = Dot @@@ Partition[tab, Divisors[3960][[43]]]; Developer`PackedArrayQ /@ slow (* {False, ...


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


16

Vectorization will help a lot: a[x_?NumericQ] := N[Exp[-Abs[x]]]; x = Table[-10 + 0.02 (j - 1), {j, 1, 1001}]; A = Outer[a[#1 - #2] &, x, x]; // AbsoluteTiming (* {2.11988, Null} *) B = Exp[-Abs[x - #]] & /@ x; // AbsoluteTiming (* {0.016182, Null} *) A == B (* True *) Notice that I am doing arithmetic on vectors the size of x instead of ...


15

It turns out ListSurfacePlot3D does a terribly poor job of approximating the surface in the OP, otherwise one will just apply DiscretizeGraphics to the output obtained from ListSurfacePlot3D and be done with it. But since that's not applicable here, we present an approach that uses alpha shapes to approximate the shape of the given point set by tuning a ...


15

Partition[ Array[a, 4], 2] will do it. In general, makeMat[n_, m_] := Partition[ Array[a, n*m], m]


14

You are needlessly computing exactly the same BesselJ and BesselJZero function values over and over again. As an example, in the 100x100 case, BesselJZero[1 + l, 1] (l is a constant) is computed 551 times! You just need to compute each one once. The easy way to do that is to memoize: bj[n_?NumericQ, z_?NumericQ] := bj[n, z] = N[BesselJ[n, z]] ...


14

Simplify[MapThread[Max, {mat, 0 mat}, 2], Assumptions -> {a > 0, b > 0, d > 0}] (* {{a, 0, 0}, {d, a, 0}, {0, 0, 0}} *) Simplify[MapThread[Max, {-mat, 0 mat}, 2], Assumptions -> {a > 0, b > 0, d > 0}] (* {{0, 0, b}, {0, 0, 0}, {0, 0, a}} *)


14

Using an undocumented function: mat /. x_?Internal`SyntacticNegativeQ :> 0 {{a, 0, 0}, {d, a, 0}, {0, 0, 0}} % - mat {{0, 0, b}, {0, 0, 0}, {0, 0, a}}


13

As Daniel Lichtblau suggested in the comment, ArrayFlatten is the way. ArrayFlatten[{{A1, A2}}] // MatrixForm gives your Out[198] ArrayFlatten[{{A1}, {Transpose @ A2}}] // MatrixForm gives your Out[201] ArrayFlatten[{{A11, A12}, {A21, A22}}] // MatrixForm gives your Out[206]


13

Outer is highly optimized for several built-in functions (Plus, Times, List). Therefore Exp@-Abs@Outer[Plus, #, -#] &@Range[-10, 10, 0.02]; // RepeatedTiming (* {0.025, Null} *) gives ~50x speedup over Outer[#1 - #2&, #, #] and ~15x speedup over Outer[Subtract, #, #]. Also is a bit faster then Kuba's Exp[-Abs[x - # & /@ x]].


13

Lots of solutions. Time for a benchmark. My own contribution is Part: m = {{{1}, {2}, {3}}, {{2}, {4}, {6}}}; m[[All, All, 1]]; {{1, 2, 3}, {2, 4, 6}} Update: I made a complete mess of my earlier attempt at benchmarking. Here is a rewrite. methods = Hold[Flatten /@ m, ArrayReshape[m, Most@Dimensions@m], ArrayReshape[m, Dimensions[m][[1 ;; ...


13

You need to understand that Mathematica prefers to write some numbers in their closed form because with numerical values, you would loose information and probably precision. It is kind of why it is better to keep Sin[4] and not use -0.756802, because Sin[4] can probably later in your calculation be combined or simplified with other expressions. That being ...


12

Now fixed in version 10.2. In[1]:= m = {{0, 1}, {-1, 0}}; In[2]:= {AntihermitianMatrixQ[m], HermitianMatrixQ[m], AntihermitianMatrixQ[m]} Out[2]= {True, False, True} As per the comments, yes, there is information stored in the internal representation of matrices (for example, a symmetry flag) and no, it is not accessible from top level code.


12

I needed this decomposition to answer another question, so I broke down and implemented it myself. The code is more or less a straightforward translation of the pseudocode in Golub/Van Loan: LDLT[mat_?SymmetricMatrixQ] := Module[{n = Length[mat], mt = mat, v, w}, Do[ If[j > 1, w = mt[[j, ;; j - 1]]; v = ...


12

The problem is that Dot does not evaluate if one of the arguments is not a list of some sort, and so when the matrix A0 is added, it does something you might not foresee. Consider: m = {{1, t}, {t, -1}}; A0 = {{1, 0}, {0, 1}}; m.A[t] + A0 (* { { 1 + {{1, t}, {t, -1}}.A[t], {{1, t}, {t, -1}}.A[t] } , { {{1, t}, {t, -1}}.A[t], 1 + {{1, t}, {t, -1}}.A[t] } ...


12

Setting the Method option to "CofactorExpansion" results in the correct output. mat = {{2, 2.161209223472559` + 1.682941969615793` I}, {2.161209223472559` - 1.682941969615793` I, 2}} Inverse[mat, Method -> "CofactorExpansion"] $\ $ {{-0.57092 + 0. I, 0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092 + 0. I}} As you want to perform ...


12

Wrap the matrix rows with the the Flatten function M = {{{1}, {2}, {3}, {4}}, {{5}, {6}, {7}, {8}}} To save time you can wrap your whole matrix using: Map[Flatten, <yourmatrix> ] Map[Flatten,{{{1}, {2}, {3}, {4}}, {{5}, {6}, {7}, {8}}}] the outermost list contains two elements (the rows). the Map function wraps these elements with the ...


12

tl;dr If these functions cannot decide, they will simply return False. A False result means that the selected equality testing method wasn't able to prove equality, but it does not mean that it was able to prove inequality. Interpret the result relative to the used SameTest option value. I will try to explain what I think is happening, though some of ...


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

It's an issue of growth of term size in this example. If you do the straight iteration then at each step the matrix elements roughly quadruple in size (because each time you multiply every element by a variable, and sum four such products per matrix entry). We confirm this below on an example of half the size. tab = Partition[#, 4] & /@ ...


11

The definition of the scalar product in your question assumes that all your kets are orthogonal unit vectors. In that case, the most natural approach would be to use the built-in Bra and Ket as follows: Ket /: Dot[Bra[x__], Ket[y__]] := Times @@ MapThread[KroneckerDelta, {{x}, {y}}] BraKet[x_, y_] := Bra[x].Ket[y] Bra[2, 4].Ket[2, 4] (* ==> 1 *) ...


11

This is almost certainly an out-of-memory crash. The underlying issue is that the OS X front-end is a 32-bit program, so has a process memory limit of around 2 GB. It is normal for an attempted allocation beyond that limit to lead to a crash. A similar size ArrayPlot example I tried on my Windows machine (where the front-end is 64-bit) used more than 3 GB ...


11

As I have previously noted, QRDecomposition[] is by default set to return the so-called "thin QR" or "economy QR" decomposition; this is often the form desired in applications, since the triangular factor does not have the unneeded zero rows. MATLAB's qr(), by contrast, returns the full QR decomposition by default, and the economy QR through an option ...


11

I would simply do Dot[Conjugate[#1], #2] & @@@ Partition[mm, 2, 1] and do away with the sequential numbers, as those are easy enough to generate when needed any way (MapIndexed?) ListPlot does not need these indices. But instead of jumping straight to the solution, let's take your code and improve it step by step. Instead of arr[[i]][[j]] you can ...


11

The Pitsianis-Van Loan algorithm turns out to be surprisingly easy to implement in Mathematica: nearestKroneckerProduct[mat_?MatrixQ, dim1_?VectorQ, dim2_?VectorQ] /; TrueQ[Dimensions[mat] == dim1 dim2] := Module[{bv, cv, sig}, {bv, sig, cv} = SingularValueDecomposition[Flatten[ Map[Composition[Flatten, Transpose], ...


11

Here's a nice one-liner: TransformationFunction[{{1, 0, 0, 0}, {0, 0, -1, 0}, {0, 1, 0, 3}, {0, 0, 0, 1}}] // InverseFunction TransformationFunction[{{1, 0, 0, 0}, {0, 0, 1, -3}, {0, -1, 0, 0}, {0, 0, 0, 1}}] Note that TransformationFunction[] is the head of the results returned by geometric *Transform functions, which take a homogeneous transformation ...


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



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