CODEWITHSUNDEEP
pythonnumpymatrix
ajknzhol
ajknzhol

Reputation: 6450

Raise diagonal matrix to the negative power 1/2

I am trying to compute the matrix which has the following equation.

S = (D^−1/2) * W * (D^−1/2)

where D is a diagonal matrix of this form:

array([[ 0.59484625,  0.        ,  0.        ,  0.        ],
       [ 0.        ,  0.58563893,  0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.58280472,  0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.58216725]])

and W:

array([[ 0.        ,  0.92311635,  0.94700586,  0.95599748],
       [ 0.92311635,  0.        ,  0.997553  ,  0.99501248],
       [ 0.94700586,  0.997553  ,  0.        ,  0.9995501 ],
       [ 0.95599748,  0.99501248,  0.9995501 ,  0.        ]])

I tried to compute D^-1/2 by using numpy function linalg.matrix_power(D,-1/2) and numpy.power(D,-1/2) and matrix_power function raises TypeError: exponent must be an integer and numpy.power function raises RuntimeWarning: divide by zero encountered in power.

How to compute negative power -1/2 for diagonal matrix. Please help.

Upvotes: 6

Views: 4076

Answers (3)

Daniel Shriki
Daniel Shriki

Reputation: 41

You can do the following: 

numpy.power(D,-1/2, where=(D!=0))

And then you will avoid getting the warning: RuntimeWarning: divide by zero encountered in power

numpy will divide every value on the matrix element-wise by it's own square root, which is not zero, so basically you won't try to divide by zero anymore.

Upvotes: -1

Ashwini Chaudhary
Ashwini Chaudhary

Reputation: 250981

If you can update D(like in your own answer) then simply update the items at its diagonal indices and then call np.dot:

>>> D[np.diag_indices(4)] = 1/ (D.diagonal()**0.5)
>>> np.dot(D, W).dot(D)
array([[ 0.        ,  0.32158153,  0.32830723,  0.33106193],
       [ 0.32158153,  0.        ,  0.34047794,  0.33923936],
       [ 0.32830723,  0.34047794,  0.        ,  0.33913717],
       [ 0.33106193,  0.33923936,  0.33913717,  0.        ]])

Or create a new zeros array and then fill its diagonal elements with 1/ (D.diagonal()**0.5):

>>> arr = np.zeros(D.shape)
>>> np.fill_diagonal(arr, 1/ (D.diagonal()**0.5))
>>> np.dot(arr, W).dot(arr)
array([[ 0.        ,  0.32158153,  0.32830723,  0.33106193],
       [ 0.32158153,  0.        ,  0.34047794,  0.33923936],
       [ 0.32830723,  0.34047794,  0.        ,  0.33913717],
       [ 0.33106193,  0.33923936,  0.33913717,  0.        ]])

Upvotes: 7

ajknzhol
ajknzhol

Reputation: 6450

I got the answer by computing thro' mathematical terms but would love to see any straight forward one liners :)

def compute_diagonal_to_negative_power():
    for i in range(4):
        for j in range(4):
            if i == j:
                element = D[i][j]
                numerator = 1
                denominator = math.sqrt(element)
                D[i][j] = numerator / denominator
    return D


diagonal_matrix = compute_diagonal_to_negative_power()

S = np.dot(diagonal_matrix, W).dot(diagonal_matrix)
print(S)

"""
[[ 0.          0.32158153  0.32830723  0.33106193]
 [ 0.32158153  0.          0.34047794  0.33923936]
 [ 0.32830723  0.34047794  0.          0.33913718]
 [ 0.33106193  0.33923936  0.33913718  0.        ]]
"""

Source: https://math.stackexchange.com/questions/340321/raising-a-square-matrix-to-a-negative-half-power

Upvotes: 0

Related Questions