Tag Info

Hot answers tagged

37

Mathematica is every bit as fast as Matlab for these types of computations. The source of the discrepancy arises from the fact that Timing keeps track of total time used by all processors when Mathematica distributes the computation across them. We can examine a fair comparison using AbsoluteTiming, which is more comparable to Matlab's tic and toc. ...


35

I guess that V9 now adds this capability: $Assumptions = { Element[A, Matrices[{m, n}]], Element[B, Matrices[{n, k}]] }; TensorReduce[ Transpose[Transpose[A].Transpose[B]] ] (* Out: B.A *)


30

You're looking for ArrayFlatten. For your example matrices, R = ArrayFlatten[ {{A, {t}\[Transpose]},{0, 1}} ] (* => {{1, 0, 0, 1}, {0, 0, 1, 1}, {0, -1, 0, 1}, {0, 0, 0, 1}} *) The construct {t}\[Transpose] is necessary for ArrayFlatten to treat t as a column matrix. Then to find $\boldsymbol{R}^{-1}$, you run Inverse[R] (* => {{1, 0, 0, ...


24

Mathematica does not support this directly. You can do things of this sort using an external package called NCAlgebra. http://math.ucsd.edu/~ncalg/ The relevant documentation may be found at http://math.ucsd.edu/~ncalg/DOWNLOAD2010/DOCUMENTATION/html/NCBIGDOCch4.html#x8-510004.4 In particular have a look at "4.4.8 NCLDUDecomposition[aMatrix, Options]" ...


22

(I've been waiting for somebody to ask this question for months... :D ) Here's the Mathematica implementation of the Frobenius companion matrix approach discussed by Jim Wilkinson in his venerable book (for completeness and complete analogy with built-in functions, I provide these three): PolynomialEigenvalues[matCof : {__?MatrixQ}] := Module[{p = ...


18

This might be as good a time as any to distill the collective wisdom of Messrs. Huber, McClure, and Toad R. M. As already mentioned, there is this quantity of great interest to people in the business of solving simultaneous linear equations, called the condition number, and conventionally denoted by the symbol $\kappa$. This is usually associated with a ...


18

This answer is almost entirely about mathematics and algorithms, not Mathematica implementation. I'm not sure whether such answers are welcome on this site; I hope I haven't offended. This is a very important problem in computational algebra, but is usually stated in a more sophisticated way. The usual way that one thinks about it is to consider the ...


18

You can always modify the matrix so that the most negative eigenvalue is also the one with the largest absolute value, and hence corresponds to the first in the list returned by EigenVectors. An upper bound for the largest absolute value of any eigenvalue is the Hilbert-Schmidt norm. So you can rescale your matrix by subtracting this norm times the unit ...


18

The phase (and length) of the eigenvectors is completely undetermined unless you specify extra conditions in addition to the eigenvalue equation. Given that you don't have any additional conditions, it's not surprising that there is no well-defined way to plot the real and imaginary parts of each eigenvector component. A simple condition that makes the ...


17

I found a way to dramatically improve the performance of this algorithm by using the undocumented function SparseArray`KrylovLinearSolve. The key advantage of this function is that it seems to be a near-analog of MATLAB's pcg, and as such accepts as a first argument either: a square matrix, or a function generating a vector of length equal to the length ...


17

acl already posted the crucial information needed to solve this conundrum (i.e., the definition of Internal`CompileValues[LinearSolve]), but wishes to delete his post since he had not interpreted it to give the complete answer. Therefore I re-post the following observation along with a summary of what it means. The input, Internal`CompileValues[]; ...


17

It's just a matter of the difficulty inherent in the numerical computation of determinants. Here's what Cleve Moler has to say about determinants and characteristic polynomials in chapter 10 of his book on numerical computing: Like the determinant itself, the characteristic polynomial is useful in theoretical considerations and hand calculations, but ...


17

This bug has been fixed in V10 mat = {{7/2 - I/2, -1 + I, 1/2 + 5 I/2}, {-1 + I, 5 + I, -1 + I}, {1/2 + 5 I/2, -1 + I, 7/2 - I/2}}; Eigensystem[mat] Gives: {{6, 3 + 3 I, 3 - 3 I}, {{1, -2, 1}, {1, 1, 1}, {-1, 0, 1}}} \begin{array}{ccc} 6 & 3+3 i & 3-3 i \\ \{1,-2,1\} & \left\{1,1,1\right\} & \{-1,0,1\} \\ \end{array}


16

Mathematica offers a pretty complete set of functionality for linear algebra, and it has improved in recent versions. For example, since version 5, Mathematica has offered the generalised Schur decomposition (also known as the QZ decomposition). This certainly wasn't available in earlier versions. It handles sparse matrices and many other wrinkles. And if ...


15

Initially, Mathematica is not designed for such abstract calculations. But, Mathematica is a powerful programming language, so that one can add such functionality easily. See the following examples in related area of differential geometry: calculations in symbolic dimensions Abstract calculations


15

Here is the method I outlined. I'll illustrate on a small example where we split matrix into top and bottom halves. In[794]:= SeedRandom[1111]; halfsize = 3; mat = RandomInteger[{-4, 4}, {2*halfsize, 10}] Out[796]= {{-3, -1, 3, -3, 3, 3, 3, 3, 4, 2}, {3, 3, -3, 0, 0, 1, -2, -4, 0, -1}, {-3, 4, 3, 0, -2, 4, 3, -2, -2, -2}, {2, 2, 4, 0, -4, 4, -1, -4, ...


15

Based upon your update, you are trying to solve the system $$\mathbf{A}\vec{x} = \vec{b}$$ for $\vec{x}$, so LinearSolve is exactly what you want. Also, it has the exact form LinearSolve[A, b] that you're asking for. Internally it uses a form of Gaussian elimination to solve such systems; this is most likely a variant of LU decomposition, but other ...


15

There is no need to use Eigensystem or Eigenvectors to find the axis of a rotation matrix. Instead, you can read the axis vector components off directly from the skew-symmetric matrix $$a \equiv R^T-R$$ In three dimensions (which is assumed in the question), applying this matrix to a vector is equivalent to applying a cross product with a vector made up ...


15

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients. n = 200; shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1]; shapeVector = Flatten @ Position[Flatten @ shape, 1]; symbolArray = Array[x, {n, n}]; symbolLaplacian = ...


14

If you need to be sure that the order is correct, there is a function Eigensystem that returns a list of both the eigenvalues and -vectors, which is in the right order. {eValues, eVectors} = Eigensystem[{{2, 0}, {0, 1}}]; eValues eVectors {2, 1} {{1, 0}, {0, 1}} It's probably worth using just for the slight off-chance that Eigenvalues and Eigenvectors ...


14

Actually what you want is HermiteDecomposition. It is the integer ring form of RowReduce; that latter, while working largely over the integers for the forward elimination, actually is an echelon form over the rational field. As for efficiency, you'll have to experiment to see if it meets your needs. If not, feel free to post or send examples so I'll have ...


14

We can construct this matrix directly as a SparseArray. This allows some classes of numerical matrices to be stored as packed arrays while being combined with symbolic or exact vectors (or vice versa), so there can be storage and run-time efficiency reasons for using a SparseArray, in addition to the obvious benefit of direct construction. On the other hand, ...


14

Here is a very simple way to do it: Table[1/i! D[M, {a, i}] /. a -> 0, {i, 0, 3}] (* ==> {{{15, 0}, {0, 2}}, {{0, 1}, {1, 0}}, {{1, 5}, {-5, 0}}, {{0, 0}, {0, 0}}} *) This works even if the entries are not polynomials. If they are, you can replace the arbitrary maximum 3 in the Table index by the degree of the polynomial: Max[Exponent[M, a]] ...


14

The keyboard commands Ctrl+Enter, Ctrl+, and Tab can be used to enter this format. You can also use the menu Insert > Table/Matrix to create a table of specified size with placeholders. See Entering Tables and Matrices. Depending on the meaning of the question, this may have some bearing:


14

The rank of a matrix is typically determined by performing a Gaussian elimination and is given by the number of non-zero rows. In your second case, the large number $5.4\times 10^{12}$, when eventually used as a pivot, gives a badly conditioned matrix (myM2 is the second matrix in your question): RowReduce[myM2] RowReduce::luc: Result for RowReduce of ...


13

By default, Mathematica assumes symbols to be complex. However the elements on the main diagonal of a Hermitian matrix are necessarily real. To force Mathematica to interpret the elements on diagonal of m to be real you could replace them by their real part, i.e. m = {{Re[n], a, b, b}, {Conjugate[a], Re[n], b, b}, {Conjugate[b], Conjugate[b], Re[c], ...


13

Use Eigenvalues[mat, Cubics -> True] Eigenvectors[mat, Cubics -> True] sometimes Quartics -> Truecan be needed. or ToRadicals @ Eigenvalues[ mat] ToRadicals @ Eigenvectors[ mat] In general one cannot find roots (of higher order) polynomials in terms of radicals. The reason that Mathematica allows this option is that in general it is ...


13

There is a lot going on with your code, as already pointed out by others. But, I would like to point out a few more things, and this was much, much to long for a comment. First, I would use SparseArray to generate H, instead of the pre-allocate, fill-in method you are using, as follows: H = SparseArray[ {l_, m_} /; m >= l :> Sqrt[pp]/(Sqrt[n + l - ...


13

This appears to be a genuine bug with exact arithmetic in Eigensystem. Here is a comparison to the same calculation with real numbers, for which I use the matrix mat//N: mat = {{7/2 - I/2, -1 + I, 1/2 + 5 I/2}, {-1 + I, 5 + I, -1 + I}, {1/2 + 5 I/2, -1 + I, 7/2 - I/2}}; {vals, vecs} = Eigensystem[mat] (* ==> {{6, 3 + 3 I, 3 - 3 I}, {{1, -2, 1}, ...


13

Looking at CompilePrint[compiledGlynnAlgorithm] there are some CopyTensor in it which aren't really needed. There's also a few CoerceTensor in there when it might be faster to just coerce the integer matrix once at the beginning. By slightly adjusting the function all CopyTensor and CoerceTensor go away giving a small increase in speed: ...



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