Hot answers tagged

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

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

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


15

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


15

Start data First let us get some images. I am going to use the MNIST dataset for clarity. (And because I experimented with similar data some time ago.) MNISTdigits = ExampleData[{"MachineLearning", "MNIST"}, "TestData"]; testImages = RandomSample[Cases[MNISTdigits, (im_ -> 0) :> im], 100] Let us convince ourselves that all images have the same ...


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

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

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

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

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


11

For symbolic matrices, at dimension 12, Det switches from a memoization version of cofactor expansion to one-step row reduction. If the input is exact then denominators can be (and I believe are) removed. The approximate coefficient case is another story entirely; "exact" polynomial division will fail due to round-off error. Here is a 12x12 example with a ...


11

Notice the $2\times2$ block in the upper left portion of your matrix. This is what is known as a Jordan block. Jordan blocks are well-known to be defective; that is, they do not have a complete set of eigenvectors. Eigenvectors[{{1, 1}, {0, 1}}] {{1, 0}, {0, 0}} Since Mathematica is unable to yield a complete eigenvector set, it pads the list of ...


11

Here is a very short solution: qf = a x^2 + b y^2 + c z^2 + 2 d x y + 2 e x z + 2 f y z; 1/2 D[qf, {{x, y, z}, 2}] (* ==> {{a, d, e}, {d, b, f}, {e, f, c}} *) This is just an application of the answer to Quick Hessian matrix and gradient calculation.



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