# Tag Info

3

Just for the purpose of illustration (the comments of Guesswhoitis. and bbgogfrey). here are some ways (I prefer Outer): Using: lst = {15, 15, 1, 14, 10, 14, 4, 8, 8, 14, 11, 5, 13, 0, 5, 4} then Outer[Subtract, lst, lst] // MatrixForm Table[lst[[i]] - lst[[j]], {i, 16}, {j, 16}] // MatrixForm Partition[Subtract @@@ Tuples[lst, 2], 16] // MatrixForm ...

3

I cannot reproduce the behavior you observe in my version of Mathematica (10.2 on Win7-64), so I assume that you might be working on an older version. It would be interesting if you could add your version and platform to your question for reference. Nevertheless, in my opinion the problem seems to be that the plotting function is attempting to evaluate ...

1

Perhaps BarChart3D using the "Grid" option for ChartLayout may be more useful for your aim, e.g. mat = {{1, -1, 0, 0}, {-1, 2, 0, -1}, {0, 0, 1, -1}, {0, -1, -1, 2}}; BarChart3D[mat, ChartLayout -> "Grid", ChartLabels -> {Range[4], Range[4]}, LabelingFunction -> (Placed[Style[#, Red, Bold], Above] &), ChartStyle -> Blue] or just ...

3

This is a somewhat high-brow way of showing the Cayley-Hamilton theorem, through the power of holomorphic functional calculus. As I mentioned in this answer, one of the standard ways to define a matrix function is through a Cauchy-like construction: $$f(\mathbf A) = \frac{1}{2\pi i} \oint_\gamma f(z)\, (z \mathbf I- \mathbf A)^{-1}\,\mathrm dz$$ where the ...

7

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

2

I write this answer with the caveat that I don't have Mathematica version 9 or later (the versions which now have this very belated function built-in), but with said caveat being offset by knowing a thing or two about the function of a matrix. ;) I'd have to agree with george's take that the docs for MatrixFunction[] could probably have explained things a ...

1

This will be much faster than using rule replacement... Clip[matrix1, {5, Infinity}, {0, Infinity}] If entries are integer, Threshold[matrix1,4] s/b faster, but will fail for non-integers...

6

Just to explain... (once upon a time I was also very new to this Mathematica syntax - and often confused). The mentioned "/." replaces all elements in an expression. In addition with "/;" you can add a condition when this replacement should be done. So the solution to your problem is something like matrix1 /. x_ /; x < 5 -> 0 (do not use upper case ...

3

There are three problems with your posted code. Your objective function is a 1x1 matrix rather than a scalar. Your constraint equation is malformed. One side is a 1x1 matrix, the other is a scalar. They should both be scalars. The variables should be in a flat list. @Willinski already gave a more natural way of expressing vectors in Mathematica, so that ...

8

Your corrected code: myVect = {x, y}; offDiag = {{0, 1}, {1, 0}}; Maximize[1/2*myVect.offDiag.myVect, myVect.{2, 2} - 1 == 0, myVect] {1/16, {x -> 1/4, y -> 1/4}}

3

Here is a slightly compacted reformulation of belisarius's answer: a = Take[mat, 10, 10]; b = Take[mat, 10, -15]; c = Take[mat, -15, 10]; d = Take[mat, -15, -15]; rr = ArrayRules[d - c.SparseArray[LinearSolve[a, b]]]; detr = Det[SparseArray[rr /. (pos_ /; VectorQ[pos, IntegerQ] -> expr_) :> (pos -> C @@ pos), ...

8

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

4

It was apparently missed by the other posters that KroneckerProduct[] is built-in: KroneckerProduct[mat, Rmat]

0

A wasteful method: mat = {{2, 4, -1, 5, -2}, {-4, -5, 3, -8, 1}, {2, -5, -4, 1, 8}, {-6, 0, 7, -3, 1}}; {m, n} = Dimensions[mat]; {lu, piv} = Most[LUDecomposition[mat]] // Quiet; lm = Drop[LowerTriangularize[lu, -1], None, -1] + IdentityMatrix[m]; um = LinearSolve[lm, mat]; lm.um == IdentityMatrix[m][[piv]].mat True

0

Here's a completely different take (from my other answer, which I think justifies a separate answer) that somewhat automates the idea of using HoldForm on the variables. It works as long as you name the variables consistently. It's still not perfect, since, you need to make the List of replacements, but once that's done, you don't need to do the cumbersome ...

4

Using as basis the great resource for core numerical algorithms below, I managed to implement a compiled linear solve which doesn't call MainEvaluate (so quite fast). I needed a linear solve for an optimization where the objective function requires inverting matrices, I was hesitating to use C++, but I preferred to stay in Mathematica. Resources ...

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

0

On the basis of $A$ being defined as posted: adef[g_List, h_] := Module[{n = Length[g]}, gs = ArrayFlatten[Table[g, {n}]]; hs = ArrayFlatten[ ReplacePart[ConstantArray[0, {n, n}], {i_, i_} :> h]]; hs - gs] To illustrate: adefex[g_List, h_] := Module[{n = Length[g]}, gs = ArrayFlatten[Table[g, {n}]]; hs = ArrayFlatten[ ...

5

Here's my relatively compact implementation of Glynn's formula, which incorporates the Gray code optimization: SetAttributes[GrayCode, Listable]; GrayCode[n_Integer] := BitXor[n, BitShiftRight[n]] permanent[mat_?MatrixQ] /; Equal @@ Dimensions[mat] := Module[{b = 2^(Length[mat] - 1)}, PadRight[{}, b, {1, -1}].(Times @@@ ...

1

I removed the fd[i_, j_], and replaced it with a pure function Piecewise[{{Max[s[[#1]] - strike1, 0], #2 == ts + 1}}] & Table[fd[i, j], {i, 1, as + 2}, {j, 1, ts + 1}] is equivament to Array[fd[#1,#2] &, {as + 2, ts + 1}] recursive[as_, ts_, ds_, strike1_] := With[{s = Table[i ds, {i, 0, as + 1}]}, Array[ Piecewise[{{Max[s[[#1]] - strike1, ...

4

Clear[Hma]; n = 5; Hma[j_Integer] = {{0.1*j, 0, 0, 0.2*j, 0}, {0, 0.1*j, 0, 0.2*j, 0}, {0, 0.1*j, 0, 0, 0.2*j}, {0.1*j, 0, 0, 0, 0.2*j}, {0, 0, 0.1*j, 0.2*j, 0}}; data1 = Flatten[Table[Thread[{j, Eigenvalues[Hma[j]]}], {j, n}], 1] // Chop; With your condition data2 = Flatten[Table[Thread[{j, Union[Eigenvalues[Hma[j]], SameTest ...

0

I am not too clear about your last condition, however the following code could yield an output CLOSE to what you have in mind: Hma[j_] = {{0.1*j, 0, 0, 0.2*j, 0}, {0, 0.1*j, 0, 0.2*j, 0}, {0, 0.1*j, 0, 0, 0.2*j}, {0.1*j, 0, 0, 0, 0.2*j}, {0, 0, 0.1*j, 0.2*j, 0}}; ListPlot[Table[Evaluate[Eigenvalues[Hma[j]]], {j, 1, 5, 1}], Frame -> True, ...

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

0

Some people (see The ubiquitous Kronecker product by Van Loan) have worked on finding two matrices $A, B$ of specified size whose tensor product $A \otimes B$ is closest (in a norm) to a given (larger) matrix $C$. That is, find $A, B$ which minimize $||C-A \otimes B||$. The algorithm is based on the SVD. There is a matlab implementation somewhere. It would ...

9

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

3

How about: av = Array[Subscript[a, ##] &, {2}]; bv = Array[Subscript[b, ##] &, {2}]; KroneckerProduct[av, bv] + KroneckerProduct[-av, bv] {{0, 0}, {0, 0}}

4

For a non-symmetric real matrix you can consider using LibraryLink to speed things up. It still won't be as fast as the Total/Tr answer, but it may be useful otherwise (call this C program SumUpperTriangle.c): #include "WolframLibrary.h" DLLEXPORT int SumUpperTriangle(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res) { /* Variable ...

5

I found this question quite interesting, so I thought I would collect the answers contributed in comments for future reference and to have the question appear as answered in search. I generated a slightly bigger matrix to play with, and minimally modified the code to render it independent of the size of the matrix. I also compared timings of each method to ...

0

(This answer is incomplete, but I don't know if I'll get time to finish it soon.) LUDecomposition nominally requires its input matrix to be square, as per the documentation. As far as I can see, the output in this case is correct except in the last row: {lu, a, b} = LUDecomposition[{{2, 4, -1, 5, -2}, {-4, -5, 3, -8, 1}, {2, -5, -4, 1, ...

2

A neater answer just uses Eigensystem: MapAt[Max, Eigensystem[m], {2, All}] // Transpose For input {{1, 2}, {3, 4}}, this returns {{1/2 (5 + Sqrt[33]), 1}, {1/2 (5 - Sqrt[33]), 1}}.

2

As long as you use Set (m =) rather than SetDelayed (m :=) the matrix will not be given new values unless you reevaluate the definition of m. SeedRandom[1]; Clear[m] m = RandomReal[{0, 1}, {2, 2}] {{0.817389, 0.11142}, {0.789526, 0.187803}} m {{0.817389, 0.11142}, {0.789526, 0.187803}} m {{0.817389, 0.11142}, {0.789526, 0.187803}} m ...

2

As an afterthought to my comment - if speed is important, this should handily beat existing answers, particularly on larger cases: Fold[If[#2[[1]] < #1[[-1, 1]], #1, Append[#1, #2]] &, {First@rawdata},Rest@rawdata] and this will be even faster: FixedPoint[Pick[#, UnitStep@Differences[Prepend[#[[All, 1]], 0]], 1] &, rawdata] finally, fastest ...

7

Dot @@ (a /@ Range@10) (* {{1, 55}, {0, 1}} *) Also Dot @@ Array[a, 10]

6

Here are some options for the example matrix mat = {{0, 0, 0, 1, 2, 3, 1, 3, 0, 2} , {1, 0, 0, 3, 2, 3, 1, 2, 0, 0} , {0, 1, 0, 1, 0, 3, 1, 2, 3, 2}}; The following returns a list at each spot with the binary digits as elements: list = Map[IntegerDigits[#, 2, 2] &, mat, {2}] (* {{{0, 0}, {0, 0}, {0, 0}, {0, 1}, {1, 0}, {1, 1}, {0, 1}, ...

3

This s/b considerably faster for large cases: #[[Union[SparseArray[#]["NonzeroPositions"][[All, 1]]]]] &@array and this is even faster: Replace[#, ConstantArray[0, Length@#[[1]]] -> Sequence[], {1}] &@array and about the same as latter: DeleteCases[#, ConstantArray[0, Length@#[[1]]]]&@array

0

In the absence of specific definitions for the matrices in the question, I simply made up my own. Here is a way you can make it work: J = {{1, 2}, {2, 1}}; T = 1; γi = {{1, 0}, {0, 0}}; g = Identity; Clear[y, Y] y = Function[t, ##] &[Array[Y[##][t] &, Dimensions[J]]]; eqn1 = Thread /@ Thread[y'[t] + (J\[Transpose] + ...

1

In five days no one gave an answer, so I will post what I developed, although it is a poor solution: You can fill zeros of the tensor with variables that are not used, like z1,z2,z3,.... Now the derivative over this variables is zero, so I got desired result

6

Also Most@MapThread[Dot, {Conjugate@mm , RotateLeft@mm}] Or Dot @@@ Transpose@{Conjugate@mm , RotateLeft@mm} // Most Or something like shooting yourself in the foot :D ListCorrelate[{Conjugate, Identity}, mm, 1, 0, #1[#2] &, Dot, 1] // Most

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

7

You could use: Select[m,#!={0,0,0}&] Where #!={0,0,0}& is a pure function that returns True for any list not equal to the list of three 0's.

4

If you call Orthogonalize at the end, you're orthogonalizing the eigenvectors in a different order (i.e. after sorting on eigenvalue, rather than before). Orthogonalizing the same list in a different order usually gives a different output. Orthogonalize[{{1., 2}, {1, 3}}] (* {{0.447214, 0.894427}, {-0.894427, 0.447214}} *) Orthogonalize[{{1, 3}, {1., 2}}] ...

1

Using ReplaceRepeated (//.) data1 = {{611.011, 1008}, {611.062, 1077}, {611.114, 1193}, {610.958, 894}, {611.009, 1621}, {611.061, -166}, {611.112, 704}, {611.164, 131}, {611.215, 1306}, {692.637, 6394}, {692.688, 6369}, {692.739, 6664}, {692.328, 6790}, {692.379, 7378}, {692.431, 5761}, {692.482, 6750}, {692.533, 6348}, {692.584, 7535}, ...

5

Not all parts of your file are recognized as numeric data during the Import. u = Import["file.dat"] Head /@ Flatten[u] $\$ {String, Integer, Integer, String} You can convert it for example with SetAttributes[stringToNumber, Listable]; stringToNumber[s_String] := ToExpression[StringReplace[s, {"e" -> " 10^", "i" -> " I"}]]; stringToNumber[s_] := ...

4

r1 = Pick[r, Thread[# >= FoldList[Max, #]]] &@r[[All, 1]] ListLinePlot@r1

6

One possibility: try this (you can use Import on your file instead): lines = ImportString[ "1+1e-18i 24 42 23.43e-23i", "Lines"] Then: result = Map[StringSplit, lines] And: Map[Interpreter["ComplexNumber"], result, {2}] This gives: {{1. + 1.*^-18 I, 24}, {42, 0. + 2.343*^-22 I}}

1

You can use list1 = {storagematrix[[i + 4]], storagematrix[[i + 5]]}; list2 = {a[[i]], a[[i + 1]]}; t=KroneckerProduct[list1, list2]; and Position[t,Min[t]] can obtain the indices of what you want to know about what product is less than others.

5

As noted in the docs for LUDecomposition[], the two triangles are by default returned together as a single array; this is customary for LU decomposition routines, as in the original LINPACK and MATLAB's lu(). In fact, exactly this same format is stored internally by the LinearSolveFunction[] returned by LinearSolve[]: a = {{1, 2, 3}, {2, 4, 1}, {2, 5, 7}}; ...

8

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

Top 50 recent answers are included