Hot answers tagged

22

The comments by both Michael E2 and J. M. ♦ are already an excellent answer, so this is just my attempt at summarizing. Undocumented means just what it says: there need not be any reference pages or usage messages, or any other kind of documentation. There are many undocumented functions and if you follow MSE regularly, you will encounter them often. Using ...


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


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

How should I know the correct order of arguments without trying several times? You can't, usually. A lot of the undocumented usage that you see on this site will have been worked out by trial and error. Sometimes it is fruitless - I have explored plenty of interesting-sounding internal functions and got nowhere. Are there detailed usage information of ...


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

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


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

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

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

You can use Reduce[] to find a set of all conditions as follows: A = {{-1, 1, -1, 0, 0, 0}, {1, 0, 0, -1, -1, 0}, {0, -1, 0, 0, 1, -1}, {0, 0, 1, 1, 0, 1}}; b = {b1, b2, b3, b4}; x = {x1, x2, x3, x4, x5, x6} allConditions=Reduce[A.x == b, x] This returns b1 == -b2 - b3 - b4 && x3 == b2 + b3 + b4 - x1 + x2 && x5 == -b2 + x1 - x4 &&...


10

As indicated in the comments, machine-precision linear algebra operations in Mathematica use the Intel MKL library optimized implementation of BLAS/LAPACK. That is the case for all platforms where MKL is available: Windows, Linux and Mac OS X (there will be no obvious MKL library files present in the layout on OS X in 10.1 or later due to static linking). ...


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

Maybe m1 = {{1, 2}, {3, 4}}; m2 = m1*10; m3 = m1/10; (plus = Inactivate[m1 + m2 + m3, Plus]) /. a_List :> MatrixForm@a Activate@plus


10

First answer (extended comment actually) You have to define better your objective function. For example, the following works: ClearAll[mat, minev] SeedRandom[1] rm = RandomInteger[10, {40, 40, 3}]; mat[t_] := N[rm.{1, t, t^2}]; minev[t_?NumericQ] := First@Eigenvalues[mat[t], -1]; Take[Table[minev[t], {t, 0, 1, .01}], 3] (* {-0.864071 - 1.30548 I, -0....


10

Mr.Wizard proposed that the slowness is due to the copying of the data. It seems that this intuition is correct. We can use the two different arguments passing mechanism in LibraryLink to test this conjecture. In LibraryLink, there are the "copied" passing and the "shared" passing. The copied passing copies the data from Mathematica to the library function ...


9

The relationship between Q and R as computed by QRDecomposition and the "full QR" results (as described by Guesswhoitis} can be found in, for instance, Wikipedia. The following illustrates how to go from the Mathematica to the Wikipedia formulation. With a as defined in the question, {q, r} = QRDecomposition[a] (* {{{1/Sqrt[5], 0, 2/Sqrt[5]}, {28/...


9

One possibility is to temporarily inactivate the arithmetic operators, like so: Block[{Times = Inactive[Times], Plus = Inactive[Plus]}, Det[{{a, b}, {-a, -b}}]] a*(-1*b)+-1*b*(-1*a)


9

It seems like RandomReal[1., 2] is automatically a packed array, whereas {1., 2.} is not. Notice; (Outer[Times, list1, Developer`ToPackedArray@{1., 2.}];) // AbsoluteTiming // First (* 0.032425 *) whereas (Outer[Times, list1, {1., 2.}];) // AbsoluteTiming // First (* 0.542086 *) Also: list1 = Developer`FromPackedArray@RandomReal[1., 1000000]; (Outer[...


9

After some investigation I can mostly explain the (correct) behavior. First off, we set this up using machine doubles. After computing the eigensystems we reorder by eigenvalues. While there is the possibility of messing up conjugate pairs (wherein right and left eigenvalues get ordered differently due to small numeric discrepancies), I checked and this ...


9

How about this furious[a_, b_] := Module[{a1, a2, a3, b1, b2, b3, c}, {a1, a2, a3} = Transpose[a, {2, 3, 4, 1}]; {b1, b2, b3} = Transpose[b, {2, 3, 4, 1}]; c = {-a3 b2 + a2 b3, a3 b1 - a1 b3, -a2 b1 + a1 b2}; Transpose[c, {4, 1, 2, 3}]] Timing results (from march's answer) for version 10.4.1 list1 = RandomReal[{-1, 1}, {32, 32, 32, 3}]; list2 = ...


8

m = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}; {lu, p, c} = LUDecomposition[m]; l = lu SparseArray[{i_, j_} /; j < i -> 1, {3, 3}] + IdentityMatrix[3]; u = lu SparseArray[{i_, j_} /; j >= i -> 1, {3, 3}]; l.u == m[[p]] (* True *) l.u is equal to a permutation of the rows of m MatrixForm /@ {l, u}


8

This is not a complete solution but should provide a start. In the special case of a square matrix we can determine conditions on rank by determining where specific numbers of eigenvalues vanish. We show how to set this up using the example from the post. mat = Normal@SparseArray[{i_, i_} -> a, 5, b] (* Out[17]= {{a, b, b, b, b}, {b, a, b, b, b}, {b, b, ...


8

In this comment, it is noted that LatticeReduce[] is now using the Nguyen-Stehle variant of LLL, so any results you might see from LatticeReduce[] can be different from a "classical" implementation of LLL. Having said this, LatticeReduce[] does take options, but through a not too transparent interface: SetSystemOptions["LatticeReduceOptions" -> ...


8

I am writing this answer only to address the Wizard's comment of "not obvious". Recall that if you take two nonzero vectors $\mathbf v_1$ and $\mathbf v_2$ as your basis, you can then represent any other vector as a linear combination of these two ($x \mathbf v_1 + y \mathbf v_2$). In Mathematica, this is equivalent to the expression Transpose[{v1, v2}].{x, ...


7

b = {{1, 1, 0, 0, 0, 0}, {1, -1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1}, {-1, -2, -3, -3, -3, -3/2}}; h = NullSpace@b (* {{0, 0, -1, 0, 1, 0}, {0, 0, -1, 1, 0, 0}} * MatrixRank@Join[b, h] (* 6 *)


7

Update: it seems this only works in 10.1 and later, but not in 10.0. This works: Assuming[(a | b | c) ∈ Matrices[{d, d}], TensorExpand@Transpose[a.(b + c)]] (* Transpose[b, {2, 1}].Transpose[a, {2, 1}] + Transpose[c, {2, 1}].Transpose[a, {2, 1}] *) We needed to tell Mathematica that a, b and c are square matrices and use TensorExpand. Check here: ...



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