Hot answers tagged linear-algebra
21
(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 = ...
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 ...
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
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
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 ...
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
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
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]]
...
11
If your matrix is diagonally dominant (in the example it is) then you can do as follows. Start with a diagonal matrix comprised of the reciprocals of the diagonal of your original matrix. Find the residual and use that to form a correction. iterate as long as needed.
Here is your example, scaled down.
ClearAll[n, s, f, aInv]
n = 1000;
s = ...
11
Here is simple (unweighted) Mma version of the Matlab implementation
of Covariance Bending.
ClearAll[covBending];
covBending[mat_, tol_: 1/10000]:=If[PositiveDefiniteMatrixQ[mat], mat,
NestWhile[(Eigensystem[#][[2]].DiagonalMatrix[
Max[#, tol] & /@(Eigensystem[#][[1]])].Transpose[
Eigensystem[#][[2]]]) &, N@mat, Min[Eigensystem[#][[1]]] < ...
11
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 ...
10
You can create your new diagonal matrix V in a single step as:
V = DiagonalMatrix@SparseArray[1/Normal[Diagonal[A]]];
On my machine, this takes 0.05 seconds, compared to 9 seconds for your code above (excluding time taken to construct A).
You can verify that they're both the same:
DiagonalMatrix[SparseArray[B]] == ...
10
There is no need for the Modulus option in CharacteristicPolynomial, since PolynomialMod serves that purpose. Assume we have a matrix m e.g. :
m = RandomInteger[10, {5, 5}]
m // MatrixForm
{{10, 1, 4, 10, 9}, {1, 9, 6, 1, 5}, {9, 7, 9, 1, 0}, {1, 10, 8, 0, 4}, {4, 0, 4, 7, 10}}
then
CharacteristicPolynomial[m, x]
2310 - 4008 x + 1739 x^2 - ...
10
As I've already mentioned in my previous answer, the method of using the QZ algorithm for matrix pencils on the Frobenius companion linearization of the polynomial eigenproblem is not always the most efficient approach. To illustrate this, I'll outline a general method for solving a hyperbolic quadratic eigenvalue problem, which is known to have all its ...
10
In years past I've taught a standard-content sophomore-level (in U.S.) linear algebra course where students used Mathematica. I know, then, that if you're interested, or if it's a requirement, you can build everything up from the simplest functions for manipulating matrices or you can directly use powerful built-in Mathematica functions (or a combination of ...
9
Not knowing anything else about the matrix, I can only suggest another alternative to the determinant (which sounds like it's prohibitively time-consuming).
If m is your matrix, try to find a root (or minimum) of
Min@Diagonal@SingularValueDecomposition[m][[2]]
instead. Unfortunately, SingularValueDecomposition is also very time consuming, but as I ...
9
Here's one quick way for polynomial matrices:
polyMat = {{15 + a^2, a + 5 a^2}, {a - 5 a^2, 2}};
Transpose[PadRight[CoefficientList[polyMat, a]], {2, 3, 1}]
{{{15, 0}, {0, 2}}, {{0, 1}, {1, 0}}, {{1, 5}, {-5, 0}}}
Alternatively (as Jens hints), you can do Flatten[PadRight[CoefficientList[polyMat, a]], {3}].
You can check that the matrix polynomial is ...
9
There might be a memory leak. More likely is that memory was used at intermediate steps and, while (I hope) freed by Mathematica, was not returned to the OS.
As for why massive memory will be consumed at all, I got into a debug version of the Mathematica kernel and had a look. By row 850 or so, in the forward elimination step, I see numbers going around in ...
9
Manipulate is probably the easiest for this specific case but here is an alternative:
DynamicModule[{select = {1, 2, 3}},
Column[{
CheckboxBar[
Dynamic[select], {1 -> "g(x)", 2 -> "y=x",
3 -> "f(x)"}],
Dynamic@Plot[Evaluate@{0.5 x + 1, 2 x - 2, x}[[select]], {x, -1, 5},
PlotRange -> {-1, 5},
AspectRatio -> 1,
...
9
Rob Knapp has some excellent publicly available notebooks for some of the things you might want. Here are three in particular that may be useful.
http://www.mathematica-journal.com/issue/v7i4/features/knapp/
http://library.wolfram.com/infocenter/Conferences/288/
http://library.wolfram.com/infocenter/Conferences/7968/
Among other things they give explicit ...
9
The determinant computation is a matter of memory use in terms of how much we want to store for subdeterminants of a Laplace expansion. Mathematica simply refuses to go that route after 11x11. YOu can do your own as below.
myDet[mat_] /; Length[mat] <= 4 := Det[mat]
myDet[mat_] :=
myDet[mat] =
Sum[mat[[1, j]]*myDet[Drop[mat, {1}, {j}]], {j, ...
9
As $P$ is explicitly constructed from eigenvectors of a self-adjoint matrix, it is unitary, i.e $P P^\dagger = I\qquad$ where the $\dagger$ is the conjugate transpose (or Hermitian conjugate, if you prefer). So, calculating the inverse is simply ConjugateTranspose[P] which is much faster than calculating it using Inverse. That said, you have to ensure that ...
8
For reference there is a built-in function CholeskyDecomposition.
For improving your existing code Array may be a minor subjective improvement:
HalfFunctionalCholesky2[matrin_List?PositiveDefiniteMatrixQ] :=
Module[{dimens, uu},
dimens = Length[matrin];
uu = ConstantArray[0, {dimens, dimens}];
Array[(uu[[#]] = makerow[matrin, #, uu, dimens]) ...
8
Generally, the reason why matrices that were supposed to be positive semi-definite but are not, is because the constraint of working in a finite precision world often introduces a wee bit of perturbation in the lowest eigenvalues of the matrix, making it either negative or complex. These errors are generally of the order of machine precision, but is enough ...
8
Having taught linear algebra using both Mathematica and Matlab, I concur with what others have said that the Mathematica's features for linear algebra include all one might need for a course in undergraduate linear algebra. Since symbolic computation is also fully integrated into Mathematica, it might be better in some ways. For example, we can solve ...
8
You were almost there. Just add the following to your code:
iC3v = Inverse /@ C3v;
sa = SolveAlways[Flatten@
Table[basis[[i]][iC3v[[k]].{x, y}] == Sum[basis[[j]][{x, y}] d[k, j, i], {j, 3}],
{i, 3}, {k, 6}],
{x, y}];
MatrixForm /@ Table[d[k, i, j], {k, 6}, {i, 3}, {j, 3}] /. sa
And you get your expected result:
$\left(
...
8
One way to do this is to use Opacity to hide a graph and empty label "" to hide a label:
Manipulate[
Plot[{0.5 x + 1, x, 2 x - 2}, {x, -1, 5},
PlotRange -> {-1, 5}, AspectRatio -> 1,
PlotStyle -> {Opacity[a], Opacity[b], Opacity[c]},
Epilog -> {
Text[If[a == 1, "f(x)", ""], {4.5, 2.7}],
Text[If[b == 1, "y=x", ""], {4.5, ...
8
For speeding up Mathematica code, a little analysis of a problem often goes a long way.
Analysis
Writing $\mathbb{G}$ for the diagonal matrix of vertex degrees and $\mathbb{A}$ for the adjacency matrix, this question seeks eigenvectors of $\mathbb{1} - \mathbb{G}^{-1/2} \mathbb{A} \mathbb{G}^{-1/2}$. That is, it looks for solutions $(\mathbf{z}, \lambda)$ ...
8
This is too long for a comment and honestly, to give a real answer, there is more information required in your question. Isn't it possible, that you give a working example, so that we see what takes long and how you implemented it?
If you are calling Eigensystem for many different input values which are know, there is still some place for speed-up. Since ...
8
A possible starting point:
With[{θ = π/3},
DeleteCases[
ParametricPlot[{{x, y}, RotationTransform[-θ][Sqrt[2] {x, y}]},
{x, -5, 5}, {y, -5, 5}, PlotRange -> {{-5, 5}, {-5, 5}}],
_Polygon, ∞]]
Only top voted, non community-wiki answers of a minimum length are eligible
