Tag Info

Hot answers tagged

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


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


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


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


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


8

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


7

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


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.


7

SparseArray[{{i_, i_} :> 2, {i_, j_} :> -1 /; Abs[i - j] == 1}, {10, 10}] // MatrixForm or SparseArray[{Band[{1, 1}] -> 2, Band[{2, 1}] -> -1, Band[{1, 2}] -> -1}, {10, 10}] // MatrixForm The following is sometimes useful: m = {{{a, b}, {b, a}}}; d = {10, ...


7

Without SparseArray: n = 10; Total[ {DiagonalMatrix[Array[-1 &, n - 1], -1], DiagonalMatrix[Array[2 &, n]], DiagonalMatrix[Array[-1 &, n - 1], 1]} ] Or strictly using Array: Array[Which[#1 == #2, 2, Abs[#1 - #2] == 1, -1, True, 0] &, {10, 10}] AND, what the heck? One more for more flexible applications: a = {2, -1, -2, 3}; n = ...


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


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


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


6

This should be faster (where n is your square dimension, e.g. 1000 for 1kX1k), easily extends to n-diagonal symmetric with no performance impact: ToeplitzMatrix[PadRight[{2, -1}, n]]


5

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


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


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


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


4

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


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


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


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


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


3

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


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


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


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


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



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