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

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

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


13

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


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]


12

To verify that matrix is a zero of its characteristic polynomial, The Characteristic polynomial of the matrix is found, then evaluated for the matrix. The result should be the zero matrix. Clear[x] a = {{-1, -4, -2}, {0, 1, 1}, {-6, -12, 2}}; n = Length[a]; p = CharacteristicPolynomial[a, x]; (Sum[ Coefficient[p, x, i] MatrixPower[a, i], {i, 0, Exponent[p, ...


12

Mathematica actually has a function purpose-built for the operation you're looking for. MatrixFunction[f, m] gives the matrix generated by the scalar function f at the matrix argument m. In your case, MatrixFunction[p, A] will return the 3-by-3 zero matrix, as desired.


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

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


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

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

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

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


11

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


11

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


10

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


10

This is really a natural fit for Outer: t = Table[{i, j}, {i, 1, 2}, {j, 1, 2}]; Outer[Apply, {Plus, Subtract, Times, Divide}, t, 2] (* ==> {{{2, 3}, {3, 4}}, {{0, -1}, {1, 0}}, {{1, 2}, {2, 4}}, {{1, 1/2}, {2, 1}}} *)


10

There is another option, using the relatively new tensor capabilities of Mathematica. This is pretty much copied from another answer by jose, but I don't need any assumptions here: TensorExpand[KroneckerProduct[X, X] + KroneckerProduct[-X, X]] (* ==> 0 *) TensorExpand[KroneckerProduct[2 X, 3 Y]] (* ==> 6 KroneckerProduct[X, Y] *) There is a ...


10

input = {{1, 1}, {1, 2}, {1, 3}, {2, 3}, {3, 3}, {4, 4}}; ConnectedComponents@Graph[UndirectedEdge @@@ input] (* {{3, 1, 2}, {4}} *)


9

For a wide variety of applications, the cost of doing a scalar product is rarely linear in the complexity of the multiplicands. Furthermore, the complexity of a product is usually larger than the complexity of the inputs. This can range from the simple case of multiplying two $n$-bit integers to get a $2n$-bit sum, to the horrible case of multiplying two ...


9

Add this to your notebook or init file $PrePrint = If[MatrixQ[#], MatrixForm[#], #] &; Then all matrices will automatically display as MatrixForm and If you want to format lists as column vectors also, try $PrePrint = Which[MatrixQ[#], MatrixForm[#], VectorQ[#], ColumnForm[#], True, #] &; Now also


9

Note: the answer below referred to a previous version of the question You may have more success with SparseArray. For instance: SparseArray[ {{i_, i_} -> 1, {i_, j_} /; Abs[i - j] == 1 -> 2}, {20000, 20000}, 0 ]; // RepeatedTiming (* Out: {0.21, Null} *) You would use patterns to assign values to positions within the matrix determined ...


9

Based on the description you provided from the MATLAB documentation, corr2 is computed as $$\frac{\sum_m \sum_n (A_{mn} - \bar{A}) (B_{mn} - \bar{B})}{\sqrt{\left(\sum_m \sum_n (A_{mn} - \bar{A})^2\right) \left(\sum_m \sum_n (B_{mn} - \bar{B})^2\right)}} $$ Assuming that the mean2 function that gives the values of $ \bar{A} $ and $ \bar{B} $ does the ...


9

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


8

TL;DR; myArrayFlatten = Flatten /@ Flatten[#, {{1, 3}}] & {x,y,z} or {{x},{y},{z}} is considered column vector. Usage example: v1 = {1, 2, 3}; vNew = {{4, 5}, {6, 7}, {8, 9}}; v4 = {10, 11, 12}; v5 = {13, 14, 15}; v6 = List /@ {16, 17, 18}; { {v1, vNew, v1}, {v4, v5, v6, v1}, {vNew, vNew} } // myArrayFlatten // MatrixForm Flatten ...


8

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} *)



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