CODEWITHSUNDEEP
matlabmatrix
0x90
0x90

Reputation: 40982

How does matlab run this code so quickly?

I have this code in matlab

b = 0.25*ones(4)
a = [0 1 1 1 ; 1/3 0 0 0 ; 1/3 0 0 0; 1/3 0 0 0] 

m = .85*a + .15*b

v = [1/4 1/4 1/4 1/4]

m^1e308*v'
  1. How does matlab run m^1e308*v' so quickly? it should mutiply the matrix 1e300 times, but it probably do some other calculation, what is it?

  2. why does m^1e309*v' gives :

    ans =

       NaN
   NaN
   NaN
   NaN

  3. how can I see what m^inf is without using symbolic variables?

Upvotes: 3

Views: 373

Answers (5)

Amro
Amro

Reputation: 124563

As others have explained, diagonalization can be used to compute the powers of a matrix efficiently. Try it yourself:

[V,D] = eig(M);             %# M = V*D*inv(V)
V*(D.^realmax)/V            %# M^n = V*D^n*inv(V)

and compare against:

M^realmax                   %# same as: mpower(M,realmax)

EDIT:

Here is the answer in MathWorks own words:

^

Matrix power. X^p is X to the power p, if p is a scalar. If p is an integer, the power is computed by repeated squaring. If the integer is negative, X is inverted first. For other values of p, the calculation involves eigenvalues and eigenvectors, such that if [V,D] = eig(X), then X^p = V*D.^p/V.

If x is a scalar and P is a matrix, x^P is x raised to the matrix power P using eigenvalues and eigenvectors. X^P, where X and P are both matrices, is an error.

Upvotes: 2

Oliver Charlesworth
Oliver Charlesworth

Reputation: 272467

How does matlab run m^1e308*v' so quickly?

I can't tell you what the internals of Matlab are doing, but note that in general, A^n can be done in O(log n) time, not O(n) time.

For example, A^16 = (((A^2)^2)^2)^2.

You can also use the eigendecomposition to turn this into scalar powering, i.e. if A = U*V*U', then the power is U * V^N * U', where V is a diagonal matrix.

why does m^1e309*v' give NaN?

Double-precision cannot represent 1e309.

how can I see what m^inf is without using symbolic variables?

Use the eigendecomposition described above. If any of the eigenvalues are smaller than 1, they will go to 0, if any of them are greater than 1, they will go to infinity.

Upvotes: 8

kevlar1818
kevlar1818

Reputation: 3125

MATLAB is so fast because it uses a Basic Linear Algebra Subprograms (BLAS) Library. The level of optimization goes right down to the CPU level, with AMD and Intel writing optimized libraries for their architectures.

BLAS is used to build the foundation of the MATLAB (MATrix LABoratory), as BLAS was used to build Linear Algebra Package (LAPACK). LAPACK provides all the low-level matrix commands in MATLAB (such as the transpose ' and matrix multiplication in m = .85*a + .15*b). It's written in FORTRAN and is open source, so what most people don't know is that MATLAB is arguably selling you convenience and a GUI.

Upvotes: 2

Viktor Mellgren
Viktor Mellgren

Reputation: 4506

I don't know the internals of matlab, but diagonalisation might be the solution (also it might parallellize the calculation internaly?)

check this chapter: Applications of diagonalisation for details

Upvotes: 0

Jerome WAGNER
Jerome WAGNER

Reputation: 22412

Does the solution converge ? After one exponent n maybe it just can see that m^n * m = m^n and so m^1e30000=m^n

It can use even an heuristic because those probability matrices (sum in each column = 1 ; each value between 0 and 1) are known to converge at infinity.

Upvotes: 1

Related Questions