UPDATE 2
The problem might have been from how we import scipy, from what I can see. I am not able to replicate the issue anymore after explicitly importing scipy.linalg instead of all of scipy. For now I will accept tch's answer. Thanks a lot!
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
UPDATE I initially accepted the answer by tch but running their code --depending on the input-- does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:
When I try the following
n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = [email protected](lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)
And then loop over the answers from k=10 to k=49:
for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)
results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)
I sometimes get meaningless results:
Can somebody shed a light on this please?
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
UPDATE I initially accepted the answer by tch but running their code does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:
When I try the following
n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = [email protected](lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)
And then loop over the answers from k=10 to k=49:
for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)
results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)
I sometimes get meaningless results:
Can somebody shed a light on this please?
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
UPDATE I initially accepted the answer by tch but running their code --depending on the input-- does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:
When I try the following
n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = [email protected](lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)
And then loop over the answers from k=10 to k=49:
for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)
results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)
I sometimes get meaningless results:
Can somebody shed a light on this please?
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
UPDATE I initially accepted the answer by tch but running their code does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:
When I try the following
n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = [email protected](lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)
And then loop over the answers from k=10 to k=49:
for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)
results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)
I sometimes get meaningless results:
Can somebody shed a light on this please?
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
I do need to numerically calculate the following forms for any $x\in\mathbb{R}^n$, possibly in python:
- $x^T M^k x$, where $M\in\mathbb{R^{n\times n}}$ is a PSD sparse matrix, $n$ can be quite large in the orders of $10^3$ or higher, and where $k$ can get large values $-$possibly up to order of hundreds. I prefer $k$ to be a real number, but it is ok if it can only be an integer, if that makes a considerable accuracy difference.
- Similarly I am interested in $x^Te^{-tM}x$ where $t$ is a non-positive real value.
For case 1 I can either:
- Use the scipy.linalg.fractional_matrix_power to calculate $M^k$ and then derive $x^TM^Kx$, or
- Use scipy.linalg.svd to find SVD of $M$ as $U\Lambda U^T$ and then finally evaluate the desired value using $x^TU\Lambda^k U^Tx$. According to Cleve Moler et. al's paper "Nineteen dubious..." (method 6) SVD shouldn't be a bad idea, but I am still not sure for my overall calculation which is $x^T MX$ this is the best way, and also people here discouraged me.
- Finally if $k$ is integer, and again based on SVD I can calculate $\|x^TM^{k/2}\|^2$.
For case 2
- Again I can use off the shelf scipy.linalg.expm
- I can do SVD for $M$ and then go with $x^TU\exp(\Lambda) U^Tx$.
- Finally since I am only interested in $x^T \exp(M) x$, and not exactly $\exp(M)$ it self, I can consider the Taylor expansion of $x^T{\rm expm}(M)x\approx \sum_{i=0}^{l} \frac{1}{i!}x^TM^ix$ for some $l$ that controls the precision, and $x^TMx$ can be calculated based on case 1.
Can anybody guide me about what is the most precise way to calculate either of these expressions, up to hopefully machine precision? Any of these methods, or they're better solutions out there? I would be happy also with references.
P.S. I have cross-posted this question on math.stackexchange.com, mathoverflow, and stackoverflow and here is suggested as the write place to post it.
UPDATE I initially accepted the answer by tch but running their code does produce unstable results. So I wanted to share the experience and see if I am doing something wrong, or there is problem in the solution they proposed that I am getting unstable results:
When I try the following
n = 3000
U,_ = np.linalg.qr(np.random.randn(n,n))
U = U.astype(np.float64)
lam = np.exp(np.abs(np.linspace(-5,5,n).astype(np.float64) + np.random.randn(n)*0.4))
print (lam)
M = [email protected](lam)@U.T.astype(np.float64)
x = np.random.randn(n).astype(np.float64)
And then loop over the answers from k=10 to k=49:
for i in range(10, 50):
_,(a_,b_) = lanczos(M,x,i)
results.append((xexpMx - lanczos_quad_form(np.exp,a_,b_,np.linalg.norm(x)))/xexpMx)
I sometimes get meaningless results:
Can somebody shed a light on this please?
Became Hot Network Question