text:
## Difference between Z-test, F-test, and T-test

A z-test is used for testing the mean of a population versus a standard, or comparing the means of two populations, with large (n ≥ 30) samples whether you know the population standard deviation or not. It is also used for testing the proportion of some characteristic versus a standard proportion, or comparing the proportions of two populations.
Example:Comparing the average engineering salaries of men versus women.
Example: Comparing the fraction defectives from 2 production lines.

A t-test is used for testing the mean of one population against a standard or comparing the means of two populations if you do not know the populations’ standard deviation and when you have a limited sample (n < 30). If you know the populations’ standard deviation, you may use a z-test.
Example:Measuring the average diameter of shafts from a certain machine when you have a small sample.

An F-test is used to compare…

View original post 90 more words

Categories: 杂记

## Understanding Harris corner detector

Let’s first go over Harris detector a little bit. For a basic idea about Harris detector, check textbooks or opencv or blogs. For its intuition, check its precursor Movarec operator, which explains why we want to maximize the variation within a window to find a corner. The key to Harris detector is the variation of intensity within a sliding window $w(x,y)$ (with displacement $u$ in the x direction and $v$ in the y direction)

$E(u,v) = \sum_{x,y} w(x,y)[I(x+u, y+v) - I(x,y)]^2$

where:

• $w(x,y)$ is the window weighting function at position $(x,y)$
• $I(x,y)$ is the intensity at (x,y)
• $I(x+u, y+v)$ is the intensity at the moved window $(x+u, y+v)$

The aim is to maximize $E(u,v)$. Applying Taylor expansion and some arithmetic operations, we can get

$E(u,v) \approx \sum_{x,y} w(x,y) (u^2I_x^2+2uvI_xI_y+v^2I_y^2) \\ = (u,v)\left( \sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2\end{bmatrix}\right) \begin{pmatrix}u\\v \end{pmatrix} \\ =(u,v)M\begin{pmatrix}u\\v \end{pmatrix}$

where $(I_x, I_y)$ is the gradient at $(x,y)$.  Matrix $M$ is known as structure tensor of a pixel.  It actually a characterization of information of all pixels within the window. Then we calculate the eigenvalues $\lambda_1$ and $\lambda_2$ of matrix $M$ to determine if window corresponds to a corner. (I)

• If $\lambda_1 \approx 0$ and $\lambda_2 \approx 0$ then this pixel $(x,y)$ has no features of interest.
• If $\lambda_1 \approx 0$ and  $\lambda_2$ has some large positive value, then an edge is found.
• If $\lambda_1$ and $\lambda_2$ have large positive values, then a corner is found.

[Note: to avoid computing the eigenvalues which is computationally expensive, in original Harris detector, the authors suggested using the following score:

to determine if a point is a corner.]

The question is WHY? Why does the matrix $M$ have the capacity to determine the cornerness? Here, we look around to the covariance matrix $\Sigma$ of gradient $(I_x, I_y)$ at a pixel. Let $\mu_{x}=E[ I_{x} ],\mu_{y}=E[ I_{y}]$. By definition of covariance matrix, we have

$\begin{array}{rl} \Sigma & = \begin{bmatrix} E[ (I_{x}-\mu_{x})(I_{x}-\mu_{x}) ] & E[ (I_{x}-\mu_{x})(I_{y}-\mu_{y}) ] \\ E[ (I_{y}-\mu_{y})(I_{x}-\mu_{x}) ] & E[ (I_{y}-\mu_{y})(I_{y}-\mu_{y}) ] \\ \end{bmatrix} \\ & = \begin{bmatrix} E[ I_{x}I_{x} ] - \mu_{x}^{2} & E[ I_{x}I_{y} ] - \mu_{x}\mu_{y}\\ E[ I_{y}I_{x} ] - \mu_{y}\mu_{x} & E[ I_{y}I_{y} ] - \mu_{y}^{2}\\ \end{bmatrix} \\ & \triangleq \begin{bmatrix} E[ I_{x}I_{x} ] & E[ I_{x}I_{y} ]\\ E[ I_{y}I_{x} ] & E[ I_{y}I_{y} ]\\ \end{bmatrix}\\ \end{array}$

The last step is derived based on the the assumption that the mean of derivatives $I_x$ and $I_y$ is zeros, i.e. $(\mu_{x},\mu_{y})=0$.

On the other hand, the structure tensor $M$ is as follows if we apply an averaging filter over the time dimension (= we give uniform weight over the pixels within the pixels):

$\begin{array}{rl} M & = \frac{1}{f}\sum_{x,y} \begin{bmatrix} I_{x}^2 & I_{x}I_{y}\\ I_{x}I_{y} & I_{y}^2 \\ \end{bmatrix} \\ & = \begin{bmatrix} \widehat{I_{x}^2} & \widehat{I_{x}I_{y}}\\ \widehat{I_{x}I_{y}} & \widehat{I_{y}^2} \\ \end{bmatrix} \\ \end{array}$

where $f$ is total number of pixels within the window, and $\widehat{f}$ is the average of $f$.

It can be seen that structure tensor $M$ is actually an unbiased estimate of the covariance matrix of the gradients for the pixels within the window. Therefore, any findings/conclusions about covariance matrix could be used to analyze structure tensors.

Let’s look at an example showing how the derivatives are distributed for different types of image patches (i.e. the window).  The following figures are adapted from Robert Collins‘s slides. In this example, the image patches are from three different types: linear edge, flat area and corner. You can image that we are testing if the pixel at the center is a corner (feature point).

The first step is calculating the derivatives $I_x$ and $I_y$ for each pixel in the image, then the distributions of ($I_x$, $I_y$) are shown in the figure below. Obviously, the distributions are discriminative fro these three different type of pixels (on the edge, on flat area and a corner).(II)

• For corners, the ranges in both $I_x$ and $I_y$ are large.
• For points on an edge, one derivative has wide distribution while the other is almost all near zeros.
• For points of a flat area, both derivatives are around zeros.

This observation is indeed the reason behind (I) we discussed above.

In Relationship between eigenvectors and directions of data distribution, we have shown that the eigenvectors of a covariance matrix are the same as the directions of principle component (i.e. the axes of a ellipsoid enclosing the data points), and the eigenvalues measure the variances of date points along eigenvectors. (Click Read more to see the relationship between eigenvalues and the semi-axes) Therefore, if the smallest eigenvalue is big enough, then it means the data have large variance along all different directions–it’s a corner!

In a word, the structure tensor is actually the covariance matrix of gradients for pixels around the pixel investigated, the the reason for current eigenvalue-based analysis is based on the property of covariance matrices.

## [OPENCV] Feature detectors and descriptors

February 16, 2013 1 comment

For a simple example of how to setup and use them, see [OPENCV] OpenCV+Eclipse|Qt|VS2010

For many computer vision tasks, the first step is usually extracting features from images. It consists of two parts: feature detection and description. By feature detection, we mean to find the pixels (or regions, objects) of interest in an image, which are represented by cv::KeyPoint. They are corresponding to the corners of objects in the image. OpenCV implements many popular feature points detection methods, such as SIFT, SURF, ORB, GoodFeatureToTrack and so on.

After these feature points are located, the next step is to describe them. An intuitive description of a feature point could be its intensity/color value. But a single intensity/color value does not provide sufficiently discriminative information so that this feature point can be uniquely matched to feature points in another image. Usually we use the information around a feature point to represent it, i.e. we use information of many pixels to represent one feature point. Based on Hubel & Wiesel’s research on information processing, orientation is the most sensitive factor to visual system. That’s why most feature techniques are based on orientation to describe feature points. For example, SIFT descriptor uses the normalized orientation histogram to describe a feature point. See figure below. Since this post is not about how the descriptors are designed, interested readers should read their papers.

To provide consistent programming interface, OpenCV designs the feature detectors and descriptor extractors following OO principle. Below is the class hierarchy of most detectors and descriptors for version 2.4.1. FeatureDetector and DescriptorExtractor are the base class for detectors and descriptors, respectively. Feature2d simplifies detection and extraction by  deriving from FeatureDetector and DescriptorExtractor.  But Feature2d is still an abstract class because it doesn’t implement the virtual function detectImpl() and computeImpl().

Classes SURF, SIFT, BRISK and ORB (and others) are concrete classes which are both detectors and descriptor extractors. In order to explicitly use them as detectors or descriptor extractors, aliases are defined. For example,

typedef SURF SurfFeatureDetector;
typedef SURF SurfDescriptorExtractor;

typedef SIFT SiftFeatureDetector;
typedef SIFT SiftDescriptorExtractor;

typedef ORB OrbFeatureDetector;
typedef ORB OrbDescriptorExtractor;

From the perspective of functionality, there’s no distinction between SurfFeatureDetector and SurfDescriptorExtractor. They are exactly the same thing except they have different names. The same is true for SiftFeatureDetecor and SiftDescriptorExtractor, and others.

From this class hierarchy diagram, it can also be seen that FAST, GFTTDetector and DenseFeatureDetector are only feature point detectors. After obtaining the feature points (std::vector<cv::KeyPoint>), we can use SURF, SIFT or any other descriptor extractors to get a description for these feature points. std::vector<cv::KeyPoint> is the data structure for data passing.

Categories: 计算机视觉, 杂记 Tags: ,

## Relationship between eigenvectors and directions of data distribution

February 7, 2013 1 comment

For approaches based on dimension reduction, such as Principle Component Analysis, one critical question is “To which direction should we project the data?” Given $N$ observations of variable $X$, we might compute the covariance matrix, and then use its eigenvectors as the direction of projection. Lots of approaches use this technique (e.g. Harris corner detector), but WHY? What’s the reason behind? This posts explains.

In the figure above, we illustrate a distribution of a 2D dataset. It’s already demeaned. The discussion below applies to any dimension size.  Suppose each data point is $N$-dimensional $X$. Intuitively, we want to find linear functions of $X$, which project $X$ into a new space. Since the data is already demeaned, $\bar{X}=0$ which means the origin of the axes is now at the center of the data,  and the linear function has no intercept part.

We start by defining the first linear function (or projection): $\alpha_1^TX$, where $\alpha_1$ is the coefficients of the linear function (see the figure above for 2D projection where $\alpha_1=(a,b)^T)$. Looking at the data plot above, the line of $\alpha_1$ should pass through the data points, and the projected points should have largest span. Mathematically speaking, we want to maximize the data variance after projecting  them to a line $\alpha$, i.e. to find $\alpha_1$ so that

$\alpha_1 = \arg\max_\alpha \{Var(\alpha^TX)\}$

Example: For 2D data, $X=(X_1, X_2)$, $\alpha=(a,b)$. $\alpha^TX=aX_1+bX_2$. It’s easy to show that the expectation $E(\alpha^TX)=aE(X_1)+bE(X_2)=a\mu_1+b\mu_2$. The variance

$\begin{array} {rl} Var(\alpha^TX) & =E[ (aX_1+bX_2)-(a\mu_1+b\mu_2) ]^2\\ & =E\left[a(X_1-\mu_1)+b(X_2-\mu_2)\right]^2 \\ & =E\left[a^2(X_1-\mu_1)^2+b^2(X_2-\mu_2)^2+2ab(X_1-\mu_1)(X_2-\mu_2)\right] \\ & =a^2Var(X_1)+b^2Var(X_2)+2abCov(X_1, X_2) \\ & =a^2\sigma_{11}+b^2\sigma_{22}+2ab\sigma_{12} \\ & =(a,b)\begin{pmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{12} & \sigma_{22} \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix}\\ & = \alpha^T\Sigma\alpha \end{array}$

where $\Sigma$ is the covariance matrix of the original data.  It is worth noticing that the finding $Var(\alpha^TX)=\alpha^T\Sigma\alpha$ can be extended to any higher dimensional variable $X$. So right now our aim becomes

$\alpha_1 = \arg\max_\alpha \alpha^T\Sigma\alpha$

Without any constraint, we could always pick a larger $\alpha$ so there’s no solution to $\alpha_1$. We choose normalized $\alpha$, which means $\|\alpha\|^2=\alpha^T\alpha=1$.

Solution: Now our problem becomes a optimization problem subject to a constraint. We use Lagrange multipliers to maximize the function

$\alpha^T\Sigma\alpha-\lambda(\alpha^T\alpha-1)$

with respect to $\alpha$. This results in

$\begin{array}{rl} \frac{d}{d\alpha}\left(\alpha\Sigma\alpha-\lambda(\alpha^T\alpha-1)\right) & =0\\ \Sigma\alpha-\lambda\alpha & = 0 \\ \Sigma\alpha & = \lambda\alpha\\ \end{array}$

Suppose the solution to this equation is $\alpha=\alpha^*$ and $\lambda=\lambda_1$.

It can be easily seen that the $\alpha_1 = \alpha^*$ is the eigenvector of $\Sigma$. When the line (represented by $\alpha$) is in the same direction as the eigenvector associated with largest eigenvalue, the variance will be the largest when we project the original data onto this line, and the variable is:

$\alpha_1^T\Sigma\alpha_1 = \alpha_1^T\lambda\alpha_1 = \lambda_1\alpha_1^T\alpha_1 = \lambda_1$

Therefore, we should choose the eigenvector of the covariance matrix with largest eigenvalue as the basis of first dimension (first principle component), and the corresponding eigenvalue is indeed the variance of that projected data in that dimension.

__________________

Following the same deduction, we can conclude that the $k$-th eigenvector $\alpha_k$ is the $k$-th principle component, and $Var(\alpha^TX)=\lambda_k$. For example, for $k=1$,

$\begin{array}{rl} &\arg\max_{\alpha_2} \alpha_2^T\Sigma\alpha_2 \\ s.t. &\alpha_2^T\alpha_2 = 1 \\ &\alpha_2\alpha_1=0 \end{array}$

(The second due to $cov(\alpha_1^TX,\alpha_2^TX)=\alpha_1^T\Sigma\alpha_2=\alpha_2^T\Sigma\alpha_1=\alpha_2^T\lambda_1\alpha_1=\lambda_1\alpha_2\alpha_1=0$.)

Use Lagrange multiplier again, you will get

$\Sigma\alpha_2 = \lambda_2\alpha_2$

___________________

If you read here, I hopeit is clear now the reason why the eigenvector with largest eigenvalue coincide with the direction of biggest variance and why the eigenvectors of a covariance matrix are used as the directions of principle components.

Note: One more question in my mind is why we are so furtunate to have such a simple yet powerful solution, but I want to stop here today.

Categories: 统计, 学习&研究 Tags: , ,

## Statistical Data Mining Basics

Statistical Data Mining Tutorials
by Andrew Moore

The following links point to a set of tutorials on many aspects of statistical data mining, including the foundations of probability, the foundations of statistical data analysis, and most of the classic machine learning and data mining algorithms.

Categories: 统计 Tags:

Categories: 杂记

## [OPENGL]Some OpenGL Examples

Here are some OpenGL programs written for course assignments, from the very basic one to some advanced application. The codes are hosted on github：https://github.com/gccheng/openGL-tutorials

1. Rotating pendulum

2. Hand & fingers

3. Solar System

4. Array operation (Rotation)

5. Quaternion algebra

6. 分段贝塞尔(Bezier)曲线


text:
## Monday, 11 August 2014

### Week of 8/4 Update

This week I've been working on on fleshing out the algorithm I described last week for computing the last invariant factor of a polynomial matrix. I've written a high-level implementation in Sage to demonstrate proof of concept. There are a few subtleties I've uncovered:
• The degree of $\bar{x}$, $N$, needs to be bigger than I originally thought. A safe bound is $N\geq{}2n*deg(A)+deg(B)+1$ although I'm not yet sure whether this bound is sharp.
• The degree of $\vec{b}$, $d_{b}$, can be fairly small, about $2\log_{q}n$. This ensures that the components are probably distinct.
• When computing $\frac{n_{i}}{d}$ such that $\frac{n_{i}}{d}\equiv{}\bar{x_{i}}(mod\;x^{N})$ there's one such for any given choice of $k=deg(d)$. However only one such choice of $k$ will actually solve the system. To figure out which choice is correct I run the extended euclidean algorithm for several rows of $\bar{x}$ simultaneously and mark values of $d$ that occur in more than one row as candidates.
• Given a candidate $d$ the other $n_{i}$'s can be determined quickly by dividing $\frac{\bar{x}_{i}}{d}$
• Any particular $\bar{x}_{i}$ has a chance of producing the wrong $d$ due to cancellation in the Cramer's rule formula $\bar{x}_{i}=\frac{|A_{i}|}{|A|}$. This occurs with probability $\frac{1}{q}$. Checking more $\bar{x}_{i}$'s for candidate $d$'s fixes this. To ensure that we get the correct $d$ as a candidate with probability $p$ we need to check $qbinom(p,n,1-\frac{1}{q})$ where qbinom is the binomial quantile function.
• Sometimes there will be no $\alpha$ such that $|A(-\alpha)|\neq{}0$, even if $A$ is non-singular (over the rational functions). In this case we can still find an answer by finding a random solution to the Block-Toeplitz system: $\begin{pmatrix}A_{d_{b}}&A_{d_{b-1}}&...&A_{0}&0&...\\A_{d_{b+1}}&A_{d_{b}}&...&A_{0}&0&...\\&&\ddots{}&&\\0&...&A_{d_{A}}&A_{d_{A}-1}&...&A_{0}\end{pmatrix}\begin{pmatrix}\bar{x}^{(0)}\\\bar{x}^{(1)}\\\vdots{}\\\bar{x}^{N}\end{pmatrix}=\vec{0}$
Next week I'll start wrapping things up as the Google Summer of Code draws to a close.

## Monday, 4 August 2014

### Week of 7/28 Update

This week I've started looking into an alternative way of computing the modulus for the polynomial ring used by SmithFormIliopolus. The problem with PolyDet is that it's painfully slow, mostly due to its heavy use of arithmetic over field extensions. The new approach would avoid the use of extensions and is based on a paper by Dixon (Exact Solutions of Linear Equations Using P-Adic Expansions (1982)). The difference is that here we're working with polynomial, not p-adic representations of integers. The goal is to compute the largest invariant factor of a matrix A. For our purposes, this is actually better than the determinant since it keeps our intermediate polynomials smaller. Here's the algorithm:
• Choose a random polynomial vector $\vec{b}$ over the ground field $\mathbb{F}$
• The degree of $\vec{b}$ (defined as the maximum of the degrees of the elements of $\vec{b}$) should be large enough to sample the range of $A$, but how big exactly remains unclear.
• Choose $\alpha$ to make $A(-\alpha)$ non-singular then let $A\leftarrow{}A(x-\alpha)$
• Solve $A\vec{x}=\vec{b}$ for the first $2n$ terms of a power series solution for $\vec{x}$ for $n$ the sum of the degrees of each row of $A$
• Perform a rational reconstruction to solve $\vec{x}=\frac{\vec{n}}{d}$ for $d$
• Take a few random linear combinations of the elements of $\vec{x}$
• For each linear combination find $f$ and $g$ of minimal degree such that $g\vec{x}\equiv{}f$ (mod $x^{2n}$)
• Take the least common multiple of the various $g$'s
• Return the polynomial $d$
Asymptotically this algorithm require about the same number of field operations as PolyDet. However, because it only works over the ground field I'm hoping it will be much faster. I'm working on a Sage example of this algorithm now to work out the details. Hopefully I'll then be able to get a Linbox implementation up and working to compare against PolyDet.

## Monday, 28 July 2014

### Week of 7/21 Update

This week I've been working on making my PolyDet code useful by creating a polynomial ring where all the operations are done modulo another polynomial. I've called this new class GivaroPolyModPoly. It reuses the Givaro Extension class which already does something similar to represent extension field elements modulo an irreducible. All I do is manually set the "irreducible" to some arbitrary polynomial. Now I can supply this ring to SmithFormIliopolus and all the intermediate polynomials will be automagically reduced modulo the polynomial I computed with PolyDet.

I've also been working on implementing parallelism in the Pascal blackbox. Some people on the FFLAS-FFPACK team have written a wrapper for task-based parallelism that we've been thinking of integrating with Linbox. The advantage to task-based parallelism is that it allows for recursive uses of parallelism, whereas things like #omp for do not. In addition the wrapper would allow us to switch to other backends such as Kaapi. I've been looking in to using these wrappers to handle parallel applyLeft() in Pascal.

## Monday, 21 July 2014

### Week of 7/14 Update

This week I focused on speeding up computations over GF(3). We have a number of matrices over GF(3) for which we want to compute the invariant factors. Fortunately, Bryan Youse has written a very efficient representation of these matrices called Sliced3. It uses 2 bits to represent each matrix element and operates over n/2 elements in parallel when computing a matrix-matrix multiplication where n is the word size. The code has been developed in parallel with the mainline Linbox source so the interface has diverged a bit from the rest of the Linbox matrix classes. Also, in order to speed up the critical applyLeft() method I had to specialize for this representation so that I use iterators directly rather than constructing temporary submatrices.

I further improved what I'm now calling my Pascal blackbox class which represents matrices of the form $A_{i,j}=c_{i+j}{i+j\choose{}j}$ with elements in GF(3). It turns out it's not necessary to be fast when computing ${n\choose{}k}$, instead you can take advantage of the fractal structure of this matrix. A $3^{n}\times{}3^{n}$ matrix $B_{i,j}={i+j\choose{}j}$ consists of 9 $3^{n-1}\times{}3^{n-1}$ submatrices, $B_{i,j}$ where $B_{1,1}=B_{1,2}=B_{1,3}=B_{2,1}=B_{3,1}$, $B_{2,2}=-A_{1,1}$ and $A_{2,3}=A_{3,3}=A_{3,2}=0$. This gives rise to a natural recursive implementation of applyLeft() which proves to be about 10 times fast when used with the Sliced3 package than my old generic sparse matrix implementation. It also only uses $O(n)$ memory instead of $O(n^{\frac{4}{3}})$

## Monday, 14 July 2014

### Week of 7/7 Update

This week I've been trying to speed up PolyDet, and to get it working with SmithFormIliopolus. Testing has revealed that, while asymptotically fast, the PolyDet code is way slower than any other code in the invariant factor pipeline. The bottleneck is at the evaluation stage, where we take $b^{2}$ polynomials of degree $d$ and evaluate them at about $bn$ points. By contrast the determinant and interpolation steps are relatively cheap. I suspect that the problem is in the GivaroPoly implementation of polynomial arithmetic. GivaroPoly and GivaroExtension extension both represent a polynomial by a std::vector. This means that a polynomial over an extension field is a vector of vectors and I suspect the doubly nested loops involved are killing performance. I'm looking into FLINT which also implements polynomial arithmetic in pure C and may be faster. Some other Linboxers have also been trying out FLINT and it sounds promising. I'm also looking into automatically switching to a Zech-log implementation which the extension field is sufficiently small, but this probably won't work for very large matrices with correspondingly large values of $d$ (roughly proportional to $\frac{n}{b}$).

I've also written a blackbox class which represents a matrix which we've been trying to get the invariant factors of. It has the structure of Pascal's triangle and the $(i,j)$ element is given by $c_{i+j}{i+j\choose{}j}$. It turns out that, with a bit of dynamic programming, ${n\choose{}k}$ can be computed in just 2 field operations using $O(n)$ space for precomputation. I took advantage of this and the regularity of the matrix to come up with a very space efficient representation that should help for the larger cases.

## Monday, 7 July 2014

### Week of 6/30 Update

This week I finished implementing the PolyDet class which computes the determinant of a matrix over a polynomial ring. The basic steps of the algorithm are:
• Given: $F=GF(p)$, $A\in{}F^{n*n}[x]$
• Let $d=\sum_{i=1}^{i=m}\max_{j\in\{1..n\}}\deg(A_{i,j})$
• (Choose $d$ to be the sum of the maximum row degrees)
• Let $d\leftarrow{}\min\left\{k|2^{k}>d+1\right\}$
• (Add 1 to $d$ and round up to the nearest power of 2)
• Let $EF=GF(p^{e})$ where $p^{e}=\min{\left\{k|\geq{}d\right\}}$
• Map the elements of $A$ to $EF$ homomorphicaly
• Choose d distinct points $p_{i}$ in $EF$
• Compute $\forall{}i\in\left\{1..d\right\}:v_{i}=det(A(p_{i}))$
• Interpolate the polynomial $P(x)$ of minimal degree that fits the points $v_{i}$
• Map the coefficients of $P(x)$ back to the ground field and return
This was a real bear to get working primarily because it involved working so much with extension fields. I used the GivaroExtension class for this and it turns out that it's still fairly buggy. A lot of algorithms in Linbox don't work for extension fields because of subtle differences in how the interfaces work. In addition, currently you can't construct an extension of an extension using GivaroExtension, you need a different extension class such as GivaroZpz which has some limitations of its own. I'm hoping to clean up some of these issues in the near future.

Having gotten polynomial determinants working the next step is to adapt the Smith Normal Form code to use this functionality. My plan is to implement a ring class that will automatically perform operations mod the determinant polynomial $P(x)$. The SmithFormIliopolus code was intended to work in a similar manner, albeit over the integers. With luck I'll be able to reuse the field extension code which already does everything mod a polynomial, and just tell it to use an arbitrary polynomial instead.

I'm also looking into using chinese remaindering to compute the invariant factors of integer matrices. With luck this should be as simple as wrapping the existing finite field code with a call to the existing CRADomain code, and performing a quick verification step. It should also provide more opportunities for parallelism since this step is embarrassingly parallel.

## Monday, 30 June 2014

### Week of 6/23 Update

This week I implemented fast polynomial interpolation and evaluation. As I mentioned last week the plan is to compute $det(A)$ for a polynomial matrix $A$ by evaluating $det(A(p))$ at a bunch of points $p$ and interpolating the polynomial that runs through those values. If the largest polynomial in $A$ has degree $d$ then we need $d+1$ points. The naive algorithm for evaluating a polynomial takes $O(d)$ time per point for a total of $O(d^{2})$ field operations. Since the polynomials from block-wiedemann have degree of order $\frac{n}{b}$ (here $n$ is the matrix dimension and very big, $b$ is the blocking factor and fairly small) we definitely needed a faster method. Although we can't do an FFT per se since our fields lack points with the requisite properties, there is a similar algorithm (Modern Computer Algebra, algorithms 10.5 and 10.11) that works for any choice of points. The new PolyInterpolation class can transform a polynomial from the coefficient domain to the point domain and back in $O(d\log^{2}{d})$ time. So far it seems to work pretty well, but it still needs a little polishing up.

I also refactored the code I'd been using to compute invariant factors and wrapped it up in its own class called CoppersmithInvariantFactors. It was getting quite complicated and this should make it significantly easier to use for other applications. I also fixed some old Linbox bugs and cleaned up a matrix representation I'd been using for testing called sparse-map-map-matrix. It's now a fully fledged SparseMatrix type.

The plan now is to automatically extend the base field as needed when doing polynomial interpolation/evaluation so that there are sufficient points to represent the original polynomial. Once I've got that working I'll pull it all together into the determinant class and I'll be able to use that for the SNF computation. My original plan was to write some kind of a wrapper for the polynomial ring GivaroPoly but now I think I might be able to use a different SNF class (SmithFormIliopoulos) that already supports doing computations modulo $det(A)$.

format: plain


text:
# author not appearing in \maketitle in amsart

I am new to the amsart document class. For some reasons \author (etc..) is not showing after I use \maketitle. This is my tex:

\documentclass[12pt,reqno,a4paper]{amsart}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{mathrsfs}
\usepackage{amsthm}
\usepackage{latexsym}
\usepackage{pictexwd}
\usepackage{pdftricks}
\usepackage{natbib}
\usepackage{setspace}
\usepackage{color}
\usepackage[affil-it]{authblk}

\begin{document}

\title{Temptation, Asymmetric information and Committment}

\author{Mark}

\affil{University of Leicester}

\date{February 2014}

\maketitle

\begin{abstract}
To be written
\end{abstract}

\section{introduction}

...

\end{document}

and this is the output:

can anybody help me out?

• Ok, I have used package \usepackage{titling} and at least I get everything but I am sure there has to be some other way...
– Matt
Jan 24, 2014 at 14:54
• Your MWE is overly complicated many irrelevant packages. Doesn't amsart have its own system for authors and their affils? Jan 24, 2014 at 14:58
• I Apologise (this is actually the first time Ihear about minimal working examples) I will do better next time.
– Matt
Jan 24, 2014 at 15:40

the way author information is handled by amsart is quite different from how it is handled by article, and the authblk package is based on the article conventions. this difference is what is causing the author information to be omitted.

if you really need the authors presented as they are in article (modified by authblk), then amsart is not going to work for you.

if you are willing to accept the amsart presentation, then you will also need to change \affil to \address or \curraddr. also, the facilities of amsthm are built into amsart, and amsmath and amsfonts are automatically loaded, so it's not necessary to specify those three packages with \usepackage.

• Hello. Yes I know about amsmath and the others. Honestly I keep them there not to forget them for when I use the article document class. Can this create any conflict? As for the author issue. I am not sure I need the same output of article. I wondered if there is a way, any way, to have the authors presented in amsart as well. Not necessarily as in article (actually I do not like the article class that's why I am trying to move to amsart...)
– Matt
Jan 24, 2014 at 15:37
• the ams document classes protect against "double loading" of these packages. however, if you want to apply some option to one of the packages, it will be ignored if added with \usepackage. Jan 24, 2014 at 15:39
• regarding author information, changing this to the style used for article would only be possible by completely redefining how the information is handled on input. that would be a major modification to the internal code of amsart.cls. there is one journal based on amsart with a more article-like topmatter block: the quarterly of applied mathematics; get the author package here. it should be easier to modify that than to start from amsart.cls itself. (please remove any ams support info if you do so.) Jan 24, 2014 at 15:53
• I see. Thank you I will look into the package and correct this as soon as I can.
– Matt
Jan 24, 2014 at 16:11


text:
# Examples¶

## Tatsu¶

The file etc/tatsu.ebnf contains a grammar for the Tatsu grammar language written in its own grammar language. It is used in the bootstrap test suite to prove that Tatsu can generate a parser to parse its own language, and the resulting parser is made the bootstrap parser every time Tatsu is stable (see tatsu/bootstrap.py for the generated parser).

Tatsu uses Tatsu to translate grammars into parsers, so it is a good example of end-to-end translation.

## Regex¶

The project examples/regexp contains a regexp-to-EBNF translator and parser generator. The project has no practical use, but it’s a complete, end-to-end example of how to implement a translator using Tatsu.

## Calc¶

The project examples/calc implements a calculator for simple expressions, and is written as a tutorial over most of the features provided by Tatsu.

## g2e¶

The project examples/g2e contains a ANTLR to Tatsu grammar translator. The project is a good example of the use of models and templates in translation. The program, g2e.py generates the Tatsu grammar on standard output, but because the model used is Tatsu‘s own, the same code can be used to directly generate a parser from an ANTLR grammar. Please take a look at the examples README to know about limitations.


text:
# American Institute of Mathematical Sciences

February  2016, 9(1): 215-234. doi: 10.3934/dcdss.2016.9.215

## On the topological characterization of near force-free magnetic fields, and the work of late-onset visually-impaired topologists

 1 Boston University, ECE Dept, 8 Saint Mary's Street, Boston, MA 02215, United States

Received  September 2014 Revised  February 2015 Published  December 2015

The Giroux correspondence and the notion of a near force-free magnetic field are used to topologically characterize near force-free magnetic fields which describe a variety of physical processes, including plasma equilibrium. As a byproduct, the topological characterization of force-free magnetic fields associated with current-carrying links, as conjectured by Crager and Kotiuga, is shown to be necessary and conditions for sufficiency are given. Along the way a paradox is exposed: The seemingly unintuitive mathematical tools, often associated to higher dimensional topology, have their origins in three dimensional contexts but in the hands of late-onset visually impaired topologists. This paradox was previously exposed in the context of algorithms for the visualization of three-dimensional magnetic fields. For this reason, the paper concludes by developing connections between mathematics and cognitive science in this specific context.
Citation: P. Robert Kotiuga. On the topological characterization of near force-free magnetic fields, and the work of late-onset visually-impaired topologists. Discrete & Continuous Dynamical Systems - S, 2016, 9 (1) : 215-234. doi: 10.3934/dcdss.2016.9.215
##### References:
 [1] F. J. Almgren and W. P. Thurston, Examples of unknotted curves which bound only surfaces of high genus within their convex hulls,, Annals of Math. $2^{nd}$ Ser., 105 (1977), 527. doi: 10.2307/1970922. Google Scholar [2] M. F. Atiyah, Mathematics in the $20^{th}$ century,, Bull. LMS, 34 (2002), 1. doi: 10.1112/S0024609301008566. Google Scholar [3] G. Burde and H. Zeischang, Knots, $2^{nd}$ ed.,, De Gruyter Studies in Math No.5, (2003). Google Scholar [4] J. C. Crager and P. R. Kotiuga, Cuts for the magnetic scalar potential in knotted geometries and force-free magnetic fields,, IEEE Trans. Mag., 38 (2002), 1309. Google Scholar [5] D. A. Ellwood, P. S. Ozsvath, A. I. Stipicz and Z. Szabo. Eds., Floer Homology, Gauge Theory, and Low-Dimensional Topology,, Clay Math. Proc. Vol.5, (2006). Google Scholar [6] J. B. Etnyre, Lectures on open book decompositions and contact structures,, in Floer Homology, (2006), 103. Google Scholar [7] G. K. Francis and B. Morin, Arnold Shapiro's eversion of the sphere,, Math. Intelligencer, 2 (1980), 200. doi: 10.1007/BF03028603. Google Scholar [8] P. Frankl and L. Pontryagin, Ein Knotensatz mit Anwendung auf die Dimensionstheorie,, Math. Annalen, 102 (1930), 785. doi: 10.1007/BF01782377. Google Scholar [9] E. Giroux, Géométrie de contact: De la dimension trois vers les dimensions supérieurs,, Proc. Int'l Congress of Mathematicians, (2002), 405. Google Scholar [10] H. Goda, Circle valued Morse theory for knots and links,, in Floer Homology, (2006), 71. Google Scholar [11] P. W. Gross and P. R. Kotiuga, Electromagnetic Theory and Computation: A Topological Approach,, MSRI Monograph No. 48, (2004). doi: 10.1017/CBO9780511756337. Google Scholar [12] R. Hiptmair, P. R. Kotiuga and S. Tordeux, Self-adjoint curl operators,, Annali di Matematica Pura ed Applicata, 191 (2012), 431. doi: 10.1007/s10231-011-0189-y. Google Scholar [13] P. R. Kotiuga, On making cuts for magnetic scalar potentials in multiply connected regions,, Jour. Appl. Phys., 61 (1987), 3916. doi: 10.1063/1.338583. Google Scholar [14] P. R. Kotiuga, An algorithm to make cuts for magnetic scalar potentials in tetrahedral meshes based on the finite element method,, IEEE Trans. Mag., 25 (1989), 4129. doi: 10.1109/INTMAG.1989.690333. Google Scholar [15] P. R. Kotiuga, Helicity functionals and metric invariance in three dimensions,, IEEE Trans. Mag., 25 (1989), 2813. doi: 10.1109/20.34293. Google Scholar [16] P. R. Kotiuga, Topology-based inequalities and inverse problems for near force-free magnetic fields,, IEEE Trans. Mag., 40 (2004), 1108. doi: 10.1109/TMAG.2004.824590. Google Scholar [17] J. C. Lagarias, J. Hass and W. P. Thurston, Area inequalities for embedded disks spanning unknotted curves,, Journ. Diff. Geom., 68 (2004), 1. Google Scholar [18] S. Smale, A classification of immersions of the two-sphere,, Trans. Amer. Math. Soc., 90 (1958), 281. Google Scholar [19] W. P. Thurston and H. E. Winkelnkemper, On the existence of contact forms,, Proc. Amer. Math. Soc., 52 (1975), 345. doi: 10.1090/S0002-9939-1975-0375366-7. Google Scholar [20] H. Whitney, Moscow 1935: Topology moving toward America,, reprinted in Hassler Whitney Collected Papers, (1992), 1. Google Scholar

show all references

##### References:
 [1] F. J. Almgren and W. P. Thurston, Examples of unknotted curves which bound only surfaces of high genus within their convex hulls,, Annals of Math. $2^{nd}$ Ser., 105 (1977), 527. doi: 10.2307/1970922. Google Scholar [2] M. F. Atiyah, Mathematics in the $20^{th}$ century,, Bull. LMS, 34 (2002), 1. doi: 10.1112/S0024609301008566. Google Scholar [3] G. Burde and H. Zeischang, Knots, $2^{nd}$ ed.,, De Gruyter Studies in Math No.5, (2003). Google Scholar [4] J. C. Crager and P. R. Kotiuga, Cuts for the magnetic scalar potential in knotted geometries and force-free magnetic fields,, IEEE Trans. Mag., 38 (2002), 1309. Google Scholar [5] D. A. Ellwood, P. S. Ozsvath, A. I. Stipicz and Z. Szabo. Eds., Floer Homology, Gauge Theory, and Low-Dimensional Topology,, Clay Math. Proc. Vol.5, (2006). Google Scholar [6] J. B. Etnyre, Lectures on open book decompositions and contact structures,, in Floer Homology, (2006), 103. Google Scholar [7] G. K. Francis and B. Morin, Arnold Shapiro's eversion of the sphere,, Math. Intelligencer, 2 (1980), 200. doi: 10.1007/BF03028603. Google Scholar [8] P. Frankl and L. Pontryagin, Ein Knotensatz mit Anwendung auf die Dimensionstheorie,, Math. Annalen, 102 (1930), 785. doi: 10.1007/BF01782377. Google Scholar [9] E. Giroux, Géométrie de contact: De la dimension trois vers les dimensions supérieurs,, Proc. Int'l Congress of Mathematicians, (2002), 405. Google Scholar [10] H. Goda, Circle valued Morse theory for knots and links,, in Floer Homology, (2006), 71. Google Scholar [11] P. W. Gross and P. R. Kotiuga, Electromagnetic Theory and Computation: A Topological Approach,, MSRI Monograph No. 48, (2004). doi: 10.1017/CBO9780511756337. Google Scholar [12] R. Hiptmair, P. R. Kotiuga and S. Tordeux, Self-adjoint curl operators,, Annali di Matematica Pura ed Applicata, 191 (2012), 431. doi: 10.1007/s10231-011-0189-y. Google Scholar [13] P. R. Kotiuga, On making cuts for magnetic scalar potentials in multiply connected regions,, Jour. Appl. Phys., 61 (1987), 3916. doi: 10.1063/1.338583. Google Scholar [14] P. R. Kotiuga, An algorithm to make cuts for magnetic scalar potentials in tetrahedral meshes based on the finite element method,, IEEE Trans. Mag., 25 (1989), 4129. doi: 10.1109/INTMAG.1989.690333. Google Scholar [15] P. R. Kotiuga, Helicity functionals and metric invariance in three dimensions,, IEEE Trans. Mag., 25 (1989), 2813. doi: 10.1109/20.34293. Google Scholar [16] P. R. Kotiuga, Topology-based inequalities and inverse problems for near force-free magnetic fields,, IEEE Trans. Mag., 40 (2004), 1108. doi: 10.1109/TMAG.2004.824590. Google Scholar [17] J. C. Lagarias, J. Hass and W. P. Thurston, Area inequalities for embedded disks spanning unknotted curves,, Journ. Diff. Geom., 68 (2004), 1. Google Scholar [18] S. Smale, A classification of immersions of the two-sphere,, Trans. Amer. Math. Soc., 90 (1958), 281. Google Scholar [19] W. P. Thurston and H. E. Winkelnkemper, On the existence of contact forms,, Proc. Amer. Math. Soc., 52 (1975), 345. doi: 10.1090/S0002-9939-1975-0375366-7. Google Scholar [20] H. Whitney, Moscow 1935: Topology moving toward America,, reprinted in Hassler Whitney Collected Papers, (1992), 1. Google Scholar
 [1] Carlos J. García-Cervera, Sookyung Joo. Reorientation of smectic a liquid crystals by magnetic fields. Discrete & Continuous Dynamical Systems - B, 2015, 20 (7) : 1983-2000. doi: 10.3934/dcdsb.2015.20.1983 [2] Serge Nicaise, Simon Stingelin, Fredi Tröltzsch. Optimal control of magnetic fields in flow measurement. Discrete & Continuous Dynamical Systems - S, 2015, 8 (3) : 579-605. doi: 10.3934/dcdss.2015.8.579 [3] Diego Rapoport. Random representations of viscous fluids and the passive magnetic fields transported on them. Conference Publications, 2001, 2001 (Special) : 327-336. doi: 10.3934/proc.2001.2001.327 [4] Naoufel Ben Abdallah, Raymond El Hajj. Diffusion and guiding center approximation for particle transport in strong magnetic fields. Kinetic & Related Models, 2008, 1 (3) : 331-354. doi: 10.3934/krm.2008.1.331 [5] Mihai Bostan. On the Boltzmann equation for charged particle beams under the effect of strong magnetic fields. Discrete & Continuous Dynamical Systems - B, 2015, 20 (2) : 339-371. doi: 10.3934/dcdsb.2015.20.339 [6] Xueke Pu. Quasineutral limit of the Euler-Poisson system under strong magnetic fields. Discrete & Continuous Dynamical Systems - S, 2016, 9 (6) : 2095-2111. doi: 10.3934/dcdss.2016086 [7] Gareth Ainsworth, Yernat M. Assylbekov. On the range of the attenuated magnetic ray transform for connections and Higgs fields. Inverse Problems & Imaging, 2015, 9 (2) : 317-335. doi: 10.3934/ipi.2015.9.317 [8] Chenxi Guo, Guillaume Bal. Reconstruction of complex-valued tensors in the Maxwell system from knowledge of internal magnetic fields. Inverse Problems & Imaging, 2014, 8 (4) : 1033-1051. doi: 10.3934/ipi.2014.8.1033 [9] Pedro Caro. On an inverse problem in electromagnetism with local data: stability and uniqueness. Inverse Problems & Imaging, 2011, 5 (2) : 297-322. doi: 10.3934/ipi.2011.5.297 [10] Frank Jochmann. A singular limit in a nonlinear problem arising in electromagnetism. Communications on Pure & Applied Analysis, 2011, 10 (2) : 541-559. doi: 10.3934/cpaa.2011.10.541 [11] Denis Serre. Non-linear electromagnetism and special relativity. Discrete & Continuous Dynamical Systems - A, 2009, 23 (1&2) : 435-454. doi: 10.3934/dcds.2009.23.435 [12] Fernando Miranda, José-Francisco Rodrigues, Lisa Santos. On a p-curl system arising in electromagnetism. Discrete & Continuous Dynamical Systems - S, 2012, 5 (3) : 605-629. doi: 10.3934/dcdss.2012.5.605 [13] Bruce Hughes. Geometric topology of stratified spaces. Electronic Research Announcements, 1996, 2: 73-81. [14] Guizhen Cui, Wenjuan Peng, Lei Tan. On the topology of wandering Julia components. Discrete & Continuous Dynamical Systems - A, 2011, 29 (3) : 929-952. doi: 10.3934/dcds.2011.29.929 [15] Fengbo Hang, Fanghua Lin. Topology of Sobolev mappings IV. Discrete & Continuous Dynamical Systems - A, 2005, 13 (5) : 1097-1124. doi: 10.3934/dcds.2005.13.1097 [16] D. Alderson, H. Chang, M. Roughan, S. Uhlig, W. Willinger. The many facets of internet topology and traffic. Networks & Heterogeneous Media, 2006, 1 (4) : 569-600. doi: 10.3934/nhm.2006.1.569 [17] Flavio Abdenur, Lorenzo J. Díaz. Pseudo-orbit shadowing in the $C^1$ topology. Discrete & Continuous Dynamical Systems - A, 2007, 17 (2) : 223-245. doi: 10.3934/dcds.2007.17.223 [18] M. Delgado-Téllez, Alberto Ibort. On the geometry and topology of singular optimal control problems and their solutions. Conference Publications, 2003, 2003 (Special) : 223-233. doi: 10.3934/proc.2003.2003.223 [19] Robert J. McCann. A glimpse into the differential topology and geometry of optimal transport. Discrete & Continuous Dynamical Systems - A, 2014, 34 (4) : 1605-1621. doi: 10.3934/dcds.2014.34.1605 [20] Juan Manuel Pastor, Silvia Santamaría, Marcos Méndez, Javier Galeano. Effects of topology on robustness in ecological bipartite networks. Networks & Heterogeneous Media, 2012, 7 (3) : 429-440. doi: 10.3934/nhm.2012.7.429

2018 Impact Factor: 0.545


text:
# Why is $\frac{d}{dx}\sin^{-1}\left(2x\sqrt{1-x^2}\right)=\frac2{\sqrt{1-x^2}}$ wrong?

In our text book's(Higher Math 1st Paper-by S U Ahmed) differentiation chapter there is a section about replacing $$x$$(inside inverse trigonometric function) with trigonometric functions. A example problem was $$\frac{d}{dx}\sin^{-1}\left(2x\sqrt{1-x^2}\right)$$ and the solution given is;

Let \begin{align*} y &=\sin^{-1}\left(2x\sqrt{1-x^2}\right)\\ &= \sin^{-1}\left(2\sin \theta \sqrt{1-\sin^2 \theta}\right)\\ &=\sin^{-1}\left(2\sin \theta \cos \theta \right)\\ &=\sin^{-1}(\sin 2\theta )\\ &=2\theta\\ &=2 \sin^{-1}x \end{align*} Now, \begin{align*} \frac{d}{dx}\sin^{-1}\left(2x\sqrt{1-x^2}\right)&=\frac{d}{dx}2 \sin^{-1}x\\ &=\dfrac{2}{\sqrt{1-x^2}} \end{align*} But plotting two functions reveals the differentiation is not actually correct. If we differentiate by parts the answer would be $$\frac{2\left(-2x^{2}+1\right)}{\sqrt{1-4x^{2}\left(1-x^{2}\right)}\sqrt{1-x^{2}}}$$

[

Now my question is why this solution is wrong?

My guess: May be this is because replacing $$x$$ with $$\sin \theta$$ changes the range of $$x$$ from $$(-\infty,\infty )$$ to $$[-1,1]$$ and may be it causes some issue.

• the range of $x$ is $[-0.7071,0.7071]$ as the domain of $\sin\alpha$ is $[-1,1]$ so $2x\sqrt{1-x^}=[-1,1]$ – Soyeb Jim Jun 13 '20 at 6:42


text:
We present a novel method for generating symbolic adversarial examples: input regions \emph{guaranteed} to only contain adversarial examples for the given neural network. These regions can generate real-world adversarial examples as they summarize trillions of adversarial examples.

We theoretically show that computing optimal symbolic adversarial examples is computationally expensive. We present a method for approximating optimal examples in a scalable manner. Our method first selectively uses adversarial attacks to generate a candidate region and then prunes this region with hyperplanes that fit points obtained via specialized sampling. It iterates until arriving at a symbolic adversarial example for which it can prove, via state-of-the-art convex relaxation techniques, that the region only contains adversarial examples. Our experimental results demonstrate that our method is practically effective: it only needs a few thousand attacks to infer symbolic summaries guaranteed to contain $\approx 10^{258}$ adversarial examples.

@misc{dimitrov2020scalable, title={Scalable Inference of Symbolic Adversarial Examples}, author={Dimitar I. Dimitrov and Gagandeep Singh and Timon Gehr and Martin Vechev}, year={2020}, eprint={2007.12133}, archivePrefix={arXiv}, primaryClass={cs.LG} }


text:
# I Are Evanescent Gravitational Waves Measurable

1. May 20, 2017

### Paul Colby

Hi,

(all discussions here are in the extreme weak field approximation about Minkowski space)

For the last couple of years I've been looking into the production and reception of radio frequency gravitational waves. It's kind of a retirement project the main goal of which is to get a better understanding of GW through futile calculations. I've gotten to the point where one of the approaches I've been looking at might actually may lead to a feasible measurement. Feasible here means likely beyond my resources but well within the existing technology.

As singularly hopeless as this may sound, it may be possible to produce and subsequently detect evanescent gravitational waves. Evanescent waves are non-propagating near fields. They occur in all wave phenomena in one guise or another. For an isolated time harmonic GW source these fields die off as $(kr)^{-5}$ whereas the radiation fields die off as $(kr)^{-1}$.

The particular example I've been considering uses garden variety quartz crystal resonators in the HF radio band. A typical 4MHz fundamental mode crystal is a quartz disk 8mm in diameter and about 0.42mm thick. There are two circular electrodes plated about 4mm in diameter plated on the front and back faces of the disk. When a sinusoidal voltage is applied a volume shear mode is excited in the plane of the disk along the "x" crystal axis (and no I don't know which way this axis points or even if it's direction is controlled by the manufacturer. Generation of gravitational waves isn't in the spec sheet).

Back to evanescent GW. The experiment would be to place two 4MHz crystals near one another. Near in this case is can-to-can which puts the crystal center to center distance at 3.45mm. A quick calculation yields,

$(kr)^{-5} = (\frac{\lambda}{2\pi r})^5 = 5\times 10^{17}$

The complete back of the envelop is,

$16 \pi G c^{-4} = 4.15\times 10^{-43} \frac{s^2}{kg\;m}$

$T_{x y} = 3.35\times 10^{6} \times Q$ Pa

This is for a drive current of 0.5 A (exceeds the spec but these things are cheep). Typical Q values run 50,000 to 100,000 so $T_{x y}$ is also a "big" number. Also I'd like to mention that $T_{x y}$ depends on the mechanical stress induced in the material and is not really a direct function of the crystal mass. The metric strain is the integral of the green function times the $T_{x y}$ so a factor on the crystal volume appears,

$V = 2\pi R^2 = 2.11\times 10^{-8} m^3$

The metric strain is approximately,

$h_{x y} = G V T k (kr)^{-5} \approx \times 10^{-22}$ [edit: I've removed Q which is in T. Q is included in T]

Working a similar estimate for the received voltage gives about 50 pico volts. This prompted me to ask people who might know good approaches to measuring at these levels. I posted this yesterday on the Electrical Engineering forum,

Signal measurement ~50 pico volts

Last edited: May 20, 2017
2. May 20, 2017

### Staff: Mentor

Aren't these near fields just the regular Newtonian gravity? Regular variations of the mass distances have been tested at various frequencies, with "0" as the easiest measurement (Cavendish).

3. May 20, 2017

### Paul Colby

No. There are no waves in Newtonian gravity evanescent or otherwise. The gravitational field of a 8mm x 0.42mm disk is quite small. The Newtonian tidal forces associated with such a disk even smaller yet. The density of quartz is 2.65 that of water if I recall correctly. That's a net mass of the disk is $5.6\times 10^{-5}$ kg.

4. May 20, 2017

### Staff: Mentor

4 MHz gives a wavelength of 75 meters, and no object moves at relativistic speeds. The Newtonian term (1/r potential) is by far the dominant part of the gravitational interaction. It doesn't matter if you call it evanescent wave or Newtonian gravity. It is there, and it is orders of magnitude stronger than everything else.

5. May 20, 2017

### Paul Colby

Could be the case. I'm not certain if in electromagnetism referring to evanescent waves as just time dependent Van der Waals forces would be very acceptable or an accurate depiction of the physics. That said, I think my shear values may be a bit unrealistic so I'll have to check my numbers.

[edit: Okay why does $k$ show up in the result and how exactly does it appear in the Newtonian calculation?]

Last edited: May 20, 2017
6. May 20, 2017

### Paul Colby

Okay, I'm going argue no, these are not time dependent tidal forces of Newtonian origin. According to the wiki

https://en.wikipedia.org/wiki/Tidal_force

The leading term goes like $r^{-3}$ and of course the speed of light doesn't appear in a Newtonian calculation so there is no opportunity for the wavelength to appear. Now there are $(kr)^{-3}$ terms along with all the others up to $(kr)^{-5}$ in the full calculation so Newtonian tidal terms are likely present but not dominant at the distances of interest to me.

7. May 20, 2017

### Staff: Mentor

I don't follow your calculations (do you have a step by step version with all the formulas before plugging in numbers?), but if you get terms that are larger than Newtonian gravity and you don't propose having a black hole handy they are wrong.

One possible cross check: What happens with your formulas if you let the frequency go to 0? Does the signal strength diverge?

8. May 20, 2017

### Paul Colby

Hum, no bias there.. What is that statement based on?

The calculation is based on the radiation integral,

$h_{\mu\nu}(t,r) = -4Gc^{-4}e^{i\omega t}\int T_{\mu\nu}(r') \frac{e^{-ik|r - r'|}}{4\pi|r-r'|} d^3r'$

which is in the Lorentz Gauge $\partial^\mu h_{\mu\nu} = 0$ since $\partial^\mu T_{\mu\nu} = 0$. This result is gauged to the transverse traceless coordinates. This is done with a 4'th order differential projector acting on the space components of $h$. That's a straight forward derivation of gauge functions $\xi_\mu$ which are harmonic (i.e. obey the scalar wave equation and have an $e^{i\omega t}$ time dependance). It's kind of long so I'll refer you to the standard text on the subject (and then ask you to tell me which one because I haven't found this anywhere for general fields as above).

There are mutually inconsistent assumptions at $\omega = 0$ and the TT-gauge projection operator I used is clearly undefined in this limit. I agree 0 would be a more comfortable (that is, the only physically acceptable) limit. Great suggestion, thanks.

9. May 20, 2017

### pervect

Staff Emeritus
I wrote a rather long post in which I estimated that the total power radiated in gravitational waves from a vibrating quartz crystal would be on the order of 3*10^-41 watts or so. This is based on an internal power flow of 10^6 watts in the crystal (I'm not positive I calculated this correctly, but the number is so low that I'd need to be off by factors of billions for it to matter.)

I won't go into more details than this unless there is some interest, because it struck me, on re-reading the thread as being off-topic. What's more to the point is this. If we assume that the far field (M/r) gives a radiated power of such a small amount, how can we possibly hope to detect terms due to the near field on the order of (M/r)^5? It would seem that (M/r)^5 should be much smaller than (M/r).

10. May 20, 2017

### Paul Colby

If $M/r$ is the order parameter then clearly not. This is tantamount to assuming $answer = 0$ in everyones mind which I question. I'm explicitly dropping the constant terms and considering only terms arising from time harmonic components. As I point out, the stress terms developed in the material don't depend on the mass. The mass appears in the sound velocity in the quartz so it's in there but not in the scale of the stress. Anyway, taking limits like $k\rightarrow 0$ have to be handled correctly. It's not clear I've done that. I need to look at the TT gauge transformation in this light. So until this is understood I don't have much more to say.

11. May 20, 2017

### Staff: Mentor

You can use the Post-Newtonian expansion, for example. For weak fields, expand in h. Alternatively, expand in terms of speed. Without black holes and without fast motion, the first order dominates.
If the signal strength diverges for $\omega \to 0$, then something is wrong. Otherwise drive your crystal with 1 µHz (doesn't even have to be resonant) and observe the bending of spacetime with the naked eye.

I think the problem is your r-5 scaling assumption. The approximation you do might be valid for the immediate vicinity of the crystal - it is not valid for the distance of 1 wavelength. Instead of scaling down to smaller length scales, we could use r-5 scaling to find the deviation at 75 meters.

$$h_{\mu\nu}(t,r) = -4Gc^{-4}e^{i\omega t}\int T_{\mu\nu}(r') \frac{e^{-ik|r - r'|}}{4\pi|r-r'|} d^3r'$$
and $T_{xy} = 3.35 MPa \cdot Q$ with $Q=100,000$ we can make a quick direct estimate. Assume that the denominator is always s=4mm, assume that the numerator is always about 1 (a very good assumption for crystal sizes much smaller than 75 m), and then just plug everything in.

$$|h_{xy}| = -4 G c^{-4} T_{xy} \frac{V}{4\pi s} \approx 6 \cdot 10^{-39}$$
WolframAlpha

That is 17 orders of magnitude smaller than your estimate in the first post, supporting my guess that the r-5 scaling was applied incorrectly.

Edit: We can also see that this contribution scales with 1/r and the strain does not depend on the frequency. While strain from tension is not directly Newtonian (not a mass as source), it is still looks very similar to Newton.
A radiative part only comes from the r'-dependent complex phase in the numerator, it is much smaller than the non-radiative term, and even that one is completely negligible.

Last edited: May 20, 2017
12. May 20, 2017

### Paul Colby

What I've done may be wrong but what you are doing I think is worse. A gauge must be chosen in which the linear strain equations,

$\rho \ddot{u}_n = \partial^mT_{n m}$ (*)

where $n$ and $m$ are summed over the spatial coordinates hold as written, or, in a gauge where these equations are modified to be consistent. This is not a small deal as the linear strain,

$S_{n m} = \partial_n u_m + \partial_m u_n$

is itself very close in form to a gauge transformation so these choices matter. Certainly the problem may be solved in any gauge but must be solved correctly in the chosen gauge. I work the TT-gauge in order that, in the absence of interatomic forces, Equations (*) or

$\ddot{u}_n=0$

holds identically. Recall I'm making the assumption,

$T_{\mu\nu}(r,t) = T_{\mu \nu}^\text{static}(r) + e^{i\omega t} T_{\mu \nu}(r)$

where I ignore entirely $T^\text{static}$. To compute the TT gauge harmonic functions, $\xi_\mu(r,t)$ must be found so that,

$h^\text{TT}_{t\mu} = h_{t\mu} - \partial_t \xi_\mu - \partial_\mu\xi_t = 0.$

Frequent use of $\partial^t \xi_\mu = ik\xi_\mu$ and the like for $\partial^t h_{t \mu} = ik h_{t \mu}$ are used. Clearly, all the spatial components must vanish as $k\rightarrow 0$. I'll slog through my derivation and try to either see how to take the limit or where I've blundered.

Last edited: May 20, 2017
13. May 21, 2017

### Paul Colby

It's pretty clear what the issue with $k\rightarrow 0$ is. The back of the envelop calculation throws the volume integral under the buss for expedience and references only the dominate radial dependence of the TT-gauge green function. Liberal use of relations like,

$h_{t,t} = \frac{1}{k^2} \partial^n\partial^m h_{nm}$

on the time harmonic components have to be used with care in the 0 frequency limit. The harmonic stress energy tensor component is also constrained by a conservation relation. Basically the bottom line is the full radiation integral with a consistent sources must be done. In the 0 frequency limit all harmonic components must vanish since the $\omega = 0$ components are already in the discarded static terms. A non-vanishing term here would clearly be a double counting.

On the one hand my not jumping to a tried and true post Newtonian formalism is justified by the extreme near field calculation I'm attempting. On the other the result must be consistent with known results. The green function I'm working with is valid to all orders in $r$ as it is a closed form expression. I've directly checked it satisfies the required scalar wave equations along with all the gauge constraints of the TT gauge. Somewhere in there is a correct answer and sure it's likely tiny.

14. May 22, 2017

Staff Emeritus
Is your mistake. Radiation is, by definition, a far field effect. This is inarguable, since it is the definition. Using a near-field approximation to estimate a far-field effect can only produce nonsense.

15. May 22, 2017

### Paul Colby

Perhaps you aren't familiar with evanescent waves? It's a common term/phenomena in many fields but I haven't seen it used in gravitation. As I've said the radiation terms are well beyond experimental reach. Here I'm raising the question about the non-radiating ones. Isn't this a fair question?

16. May 22, 2017

### Paul Colby

17. May 22, 2017

### Paul Colby

This is a really good observation. A similar issue arrises with the $k\rightarrow 0$ limit in EM, at least with the way I did it. Consider the E-field produced by 1 amp current being driven into a 100pF capacitor at 4 MHz. Now keep the current constant by arbitrarily raising the applied voltage as the frequency is reduced. In the $k\rightarrow 0$ limit this voltage diverges as does the accumulated charge on the plates. This is the EM analog of the calculation that confused me.

Also it's clear from the above example that I need to provide a better answer than I gave in post #3. Part of the field is clearly just a time dependent Newtonian one. Part is provided by GR. The question is what percentage is this GR part versus the "Cavendish" part. A number can and should be put on this value.

Last edited: May 22, 2017
18. May 30, 2017

### Paul Colby

The motivation to use the TT-gauge is to preserve the Newtonian form of continuum mechanics. The time harmonic radiation integral,

$\bar{h}_{\mu\nu}(r) = -16\pi G \int T_{\mu\nu}(r')\frac{e^{-ik|r-r'|}}{4\pi|r-r'|}d^3r'$

yields a time harmonic metric strain which is not transverse or traceless. The $T_{\mu\nu}$ in the above contains just the time harmonic part of the stress energy. I've worked through 3 methods to obtain the TT-gauge strain from the one given above. These are 1) solve for the gauge transforms 2) use the gauge invariance of the curvature which leads to,

$h^{TT}_{\mu\nu}(r) = \frac{2}{k^2}R_{\mu 4\nu 4}(r)$

which by far is the most elegant. Or, 3) my way, guess and check, arguably the least defendable. All these lead to the same 4th order differential operator. All lead to the same $(kr)^{-5}$ near field dependence.

For completeness the differential operator, $\Delta_{nmpq}$ is given where lower case latin letters are as usual restricted to 1,2,3. I use two definitions to help reduce writing,

$\Delta^2_{nm} = \epsilon_{nmp}\partial^p$

$\Delta^1_{nm} = \partial_n\partial_m+k^2\delta_{nm}$

where $\epsilon_{nmp}$ is the fully antisymmetric symbol on 3 dimensions. The full operator is,

$\Delta_{nmpq} = \frac{1}{2k^4}(\Delta^1_{np}\Delta^1_{mq}+k^2\Delta^2_{np}\Delta^2_{mq})$

The TT-gauge result is computed,

$h^{TT}_{nm}(r) = \Delta_{nmpq} (\bar{h}_{pq}(r)).$

Since for space coordinates $\eta_{nm} = \delta_{nm}$ I've played fast and loose and written everything as subscripts. Doing this with the time components can and will lead to errors.

The above completely defines where I'm getting my green function radial dependencies from. The next post will hopefully address the question in the second half of post #17.

Last edited: May 30, 2017
19. May 30, 2017

### Paul Colby

The initial motivation for this thread was a realization that the gravitational field equations appear to say something quite counter to what I had always assumed as the case. Part of the realization comes from learning the field equation in more detail and part stems from my learning some continuum mechanics, something I always considered too icky to look at.

The approach to the question first raised in #2 is to construct a mechanical example which 1) produces a time harmonic metric strain 2) produces identically 0 time harmonic strain or force in the Newtonian theory. The case initially given fails to do this because while the material density remains fixed the crystal boundary changes with time leading to a truly negligible time harmonic near field force.

Consider instead a toroidal disk of inner radius, $R_1$, and outer radius, $R_2$, fashioned from an isotropic material like steel of thickness, $\tau$. Let $u(r)e^{i\omega t}$ be the displacement from the reference configuration for a harmonic vibrational mode of angular frequency, $\omega$. The mode of interest is written in cylindrical coordinates, ${r,\phi,z}$, as

$u_{\phi}(r) = U(r)\hat{\phi}$

where $U(r)$ is the sum of two cylindrical Bessel functions such that the boundary conditions are met at $R_1$ and $R_2$ the inner and out disk radius, respectively. This mode is found by solving the Lame equation,

$\mu\nabla^2 u +(\mu+\lambda)\nabla(\nabla\cdot u) = 0$

using the separation of variables. Since the material motion is always tangent to the disk boundary, the boundary of the disk remains constant with time. Also, since $\nabla\cdot u=0$, the density of the material is constant as is the custom with shear modes. As a result there is identically 0 time dependance of the Newtonian gravitational field in and about the disk. However the gravitational field equations appear to yield a non-zero answer. This mode yields a time dependent stress component,

$T_{\phi r} = \frac{du_\phi}{dr}$

which leads to a non-vanishing metric strain with the same $(rk)^{-5}$ dependence.

20. May 30, 2017

### Staff: Mentor

I don't think physicsforums is the right place for developing the calculations. If you think it can be measured, why don't you publish it?


text:
# What does it mean if a sequence is indexed beyond its bounds?

I'm looking at a paper (On Base and Turyn Sequences by C. Koukouvinos, S. Kounias and K. Sotirakoglou) that describes an algorithm for finding specific sequences. Part of the algorithm involves finding an intermediary sequence $k$ of length $m$, $m < n$, where $n + 1$ is the length of the final sequence, and summing elements $k_1$ and $k_{n+1}$ then checking congruence. The problem is that $k$ is shorter than $n + 1$ (in some cases, much shorter). I don't think it's a typo in the paper ((33) on page 831 if you're curious), and from experimenting I don't think $n+1$ % $m$ is correct. Does anyone know what this means?

Thanks!

-
Well, I am not sure I am looking at the correct page, but from what I see, $k$ seems to carry a double index everywhere, e.g., $k_{im}$. Added: Is it possible for you to also mention the page number (between 825-837) given in the document? –  Srivatsan Aug 31 '11 at 14:03
It's on page 831 - I've edited my post to add this. The m in the double index indicates which of the intermediary sequences this is - I left it out of my description to simplify things. The n+1 only appears in the first index, so the second index didn't seem important. –  Joel Aug 31 '11 at 14:08
Well, there's no standard meaning to indexing beyond the bound, if that's your question. It could mean something in this paper, but I presume the author would've indicated it somewhere. I am not brave enough to dig deeper :-) –  Srivatsan Aug 31 '11 at 14:12
That's okay, I don't blame you :-) As far as I can tell, they don't explain it in the paper. I was hoping there was some common convention for this... –  Joel Aug 31 '11 at 14:15


text:
# Argonne National Laboratory

Seminar | LANS Seminar | Mathematics and Computer Science Division

# An Optimization-Based Perturbed Newton Method for Viscoplastic Fluids with von Mises Yielding

LANS Informal Seminar

Abstract: We target highly nonlinear fluid models with von Mises yield critera that give rise to optimization problems with Hessians exhibiting a problematic (near) null space upon linearization with Newton’s method.

The null space is caused by a projector-type coefficient in the Hessian, which is created by terms in the objective functional that resemble the $L^1$-norm. This occurs, for example, in nonlinear incompressible Stokes flow in Earth’s mantle with plastic yielding rheology, which effectively limits stresses in the mantle by weakening the viscosity depending on the strain rate. Using a standard Newton linearization for such an application is known to produce severe Newton step-length reductions due to backtracking line search and stagnating nonlinear convergence. Additionally, these effects become increasingly prevalent as the mesh is refined.

We analyze issues with the standard Newton linearization in an abstract setting and propose an improved linearization, which can be applied straightforwardly to Stokes flow with yielding and other applications as total variation regularization.

Finally, numerical experiments compare the standard and improved Newton linearizations in practice. When we employ our improved linearization within our inexact Newton-Krylov method, a fast and highly robust nonlinear solver is attained that exhibits mesh-independent convergence and scales to large numbers of cores with high parallel efficiency.


text:
# zbMATH — the first resource for mathematics

On the dimension of the zero or infinity tending sets for linear differential equations. (English) Zbl 0484.34003

##### MSC:
 34A30 Linear ordinary differential equations and systems 34D05 Asymptotic properties of solutions to ordinary differential equations
##### Keywords:
complex-valued continuous functions
Full Text:
##### References:
 [1] W. A. Coppel, Stability and asymptotic behavior of differential equations, D. C. Heath and Co., Boston, Mass., 1965. · Zbl 0154.09301 [2] Philip Hartman, The existence of large or small solutions of linear differential equations, Duke Math. J. 28 (1961), 421 – 429. · Zbl 0102.30301 [3] Philip Hartman, Ordinary differential equations, S. M. Hartman, Baltimore, Md., 1973. Corrected reprint. · Zbl 0281.34001 [4] Jack W. Macki and James S. Muldowney, The asymptotic behaviour of solutions to linear systems of ordinary differential equations, Pacific J. Math. 33 (1970), 693 – 706. · Zbl 0183.36002 [5] Marvin Marcus and Henryk Minc, A survey of matrix theory and matrix inequalities, Allyn and Bacon, Inc., Boston, Mass., 1964. · Zbl 0247.15002 [6] H. Milloux, Sur l’équation différentielle $$x'' + A(t)x = 0$$, Prace Mat. Fiz. 41 (1934), 39-53. · Zbl 0009.16402 [7] Binyamin Schwarz, Totally positive differential systems, Pacific J. Math. 32 (1970), 203 – 229. · Zbl 0193.04501 [8] M. Ō. Tnūthail, Algēbar Iolscoile, Oifig an tSolāthair, Baile Ātha Cliath, 1947.
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.


text:
## Wednesday, 4 May 2011

### Putting comments next to paragraphs in latex and emacs (fix fill-paragraph)

Again, this is just to save this for the edification of my future self.  For the rest of you, do avoid using 30 year old technology if you value achieving anything on a Wednesday.

So, I like to put little comments by my paragraphs in papers I'm writing:

% this is rubbish
blah blah blah

% this is just plain wrong
blah blah blah

% I have checked this and it is ok

Which is very useful, because this is the way to annotate the latex source so that I can understand why I've written what I did.  Now, the annoying thing is that I also like to bounce on M-q when I can't think of anything else to say, so my lovely paragraphs get turned into comments by the default emacs fill-paragraph command:

% I have checked this and it is ok debedobedo

At which point the text disappears from my output, and I get cranky.

The solution is, of course, very simple but impossible to actually find (get it together google, you are now officially beyond being useless):


text:
# Methods for Adaptive Shrinkage, using Empirical Bayes

The R package 'ashr' implements an Empirical Bayes approach for large-scale hypothesis testing and false discovery rate (FDR) estimation based on the methods proposed in M. Stephens, 2016, "False discovery rates: a new deal", . These methods can be applied whenever two sets of summary statistics---estimated effects and standard errors---are available, just as 'qvalue' can be applied to previously computed p-values. Two main interfaces are provided: ash(), which is more user-friendly; and ash.workhorse(), which has more options and is geared toward advanced users. The ash() and ash.workhorse() also provides a flexible modeling interface that can accommodate a variety of likelihoods (e.g., normal, Poisson) and mixture priors (e.g., uniform, normal).

This repository contains an R package for performing "Adaptive Shrinkage."

To install the ashr package first you need to install devtools:

To use the interior-point solver (optmethod = "mixIP"), you need to install MOSEK and the Rmosek package. We have provided some Mac-specific and Linux-specific instructions for installing MOSEK.

The main function in the ashr package is ash. To get minimal help:

## More background

The ashr ("Adaptive SHrinkage") package aims to provide simple, generic, and flexible methods to derive "shrinkage-based" estimates and credible intervals for unknown quantities $\beta=(\beta_1,\dots,\beta_J)$, given only estimates of those quantities ($\hat\beta=(\hat\beta_1,\dots, \hat\beta_J)$) and their corresponding estimated standard errors ($s=(s_1,\dots,s_J)$).

The "adaptive" nature of the shrinkage is two-fold. First, the appropriate amount of shrinkage is determined from the data, rather than being pre-specified. Second, the amount of shrinkage undergone by each $\hat\beta_j$ will depend on the standard error $s_j$: measurements with high standard error will undergo more shrinkage than measurements with low standard error.

### Methods Outline

The methods are based on treating the vectors $\hat\beta$ and $s$ as "observed data", and then performing inference for $\beta$ from these observed data, using a standard hierarchical modelling framework to combine information across $j=1,\dots,J$.

Specifically, we assume that the true $\beta_j$ values are independent and identically distributed from some unimodal distribution $g$. By default we assume $g$ is unimodal about zero and symmetric. You can specify or estimate a different mode using the mode parameter. You can allow for asymmetric $g$ by specifying mixcompdist="halfuniform".

Then, we assume that the observations $\hat\beta_j \sim N(\beta_j,s_j)$, or alternatively the normal assumption can be replaced by a $t$ distribution by specifying df, the number of degrees of freedom used to estimate $s_j$. Actually this is important: do be sure to specify df if you can.

# ashr v2.0

## Major Changes

There are a few major changes in output and input that will likely break existing dependencies. Here are the highlights

• the main output (lfsr, lfdr, etc) is rearranged into a dataframe, called result.

• so, for example, the lfsr is now a$result$lfsr instead of a\$lfsr

• or, better, use the accessor function get_lfsr(a) to extract the lfsr etc

• I added accessor functions get_lfsr, get_lfdr, get_pm, get_psd etc to access the lfsr, lfdr, posterior mean and posterior standard deviation. Using these functions to access results will help ensure your code remains valid if I happen to change the internal structure of the results again (although not anticipated...)

• output fitted.g is renamed fitted_g, and flash.data becomes flash_data to make the whole package convention more consistent. Also fit becomes fit_details.

• function prefixes comppost and compdens replaced with comp_post and comp_dens, again for consistency.

• nonzeromode option is replaced with the option mode to specify mode. Or use mode="estimate" to estimate the mode.

• more flexible control of output. For example, you can say you want only the logLR output by specifying outputlevel = c("lfsr","logLR"), or only posterior mean by outputlevel = c("PosteriorMean").

# Reference manual

install.packages("ashr")

2.2-47 by Peter Carbonetto, a year ago

https://github.com/stephens999/ashr

Report a bug at https://github.com/stephens999/ashr/issues

Browse source code at https://github.com/cran/ashr

Authors: Matthew Stephens [aut] , Peter Carbonetto [aut, cre] , Chaoxing Dai [ctb] , David Gerard [aut] , Mengyin Lu [aut] , Lei Sun [aut] , Jason Willwerscheid [aut] , Nan Xiao [aut] , Mazon Zeng [ctb]

Documentation:   PDF Manual

Imports Matrix, stats, graphics, Rcpp, truncnorm, mixsqp, SQUAREM, etrunct, invgamma

Suggests testthat, knitr, rmarkdown, ggplot2, REBayes


text:
Linear inverted pendulum model

The linear inverted pendulum model focuses on the translational components of a legged robot's dynamics.

Derivation

Both fixed and mobile robots are usually modeled as rigid bodies connected by actuated joints. The equation of motion for such a system are:

\begin{equation*} \bfM(\bfq) \qdd + \qd^\top \bfC(\bfq) \qd + \bfg(\bfq) = \bfS^\top \bftau + \sum_{i=1}^N \bfJ_{C_i}^\top \bff_i, \end{equation*}

where $$\bfq$$ is the vector of actuated and unactuated coordinates. Actuated coordinates correspond to joint angles directly controlled by motors. Unactuated coordinates correspond to the six degrees of freedom for the position and orientation of the floating base (a frame attached to any of the robot's bodies) with respect to the inertial fame. The vector $$\bfq$$ is typically high-dimensional.

Centroidal dynamics

The first working assumption to simplify this model is (Assumption 1) that the robot has enough joint torques to realize the actuated part of the equation, and focus on the Newton-Euler equations that correspond to the six unactuated coordinates:

\begin{equation*} \begin{bmatrix} m \bfpdd_G \\ \dot{\bfL}_G \end{bmatrix} = \begin{bmatrix} \bff \\ \bftau_G \end{bmatrix} + \begin{bmatrix} m \bfg \\ \boldsymbol{0} \end{bmatrix} \end{equation*}

where on the left-hand side $$\bfp_G$$ is the position of the center of mass (CoM) and $$\bfL_G$$ is the net angular momentum around the CoM, while on the right-hand side $$\bff$$ is the resultant of contact forces, $$\bftau_G$$ is the moment of contact forces around the CoM, $$m$$ is the robot mass and $$\bfg$$ is the gravity vector. This model is called centroidal dynamics.

Linearized dynamics

Angular momentum or height variations make centroidal dynamics nonlinear. This means for instance that, to generate a trajectory for this system, one needs to solve a nonlinear optimization. An alternative to linearize this system is to make two assumptions:

• Assumption 2: there is no angular momentum around the center of mass $$(\dot{\bfL}_G=\boldsymbol{0})$$. This is why the Honda P2 walks with locked arms.
• Assumption 3: the center of mass keeps a constant height. This is why the Honda P2 walks with bent knees.

Under these two assumptions, the equations of motion of the walking biped are reduced to a linear model, the linear inverted pendulum:

\begin{equation*} \bfpdd_G = \omega^2 (\bfp_G - \bfp_Z) \end{equation*}

where $$\omega^2 = g / h$$, $$g$$ is the gravity constant, $$h$$ is the constant height of the center of mass and $$\bfp_Z$$ is the position of the zero-tilting moment point (ZMP). The constant $$\omega$$ is called natural frequency of the linear inverted pendulum.

In this model, the robot can be seen as a point-mass concentrated at $$G$$ resting on a mass-less leg in contact with the ground at $$Z$$. Intuitively, the ZMP is the point where the robot applies its weight. As a consequence, this point needs to lie inside the contact surface $$\cal S$$.

To go further

The linear inverted pendulum mode (not model) was introduced in Kajita et al. (2001) for walking pattern generation. The rationale for calling it a mode of motion is that it relies solely on planar linear momentum, which is one among many ways to affect the acceleration of the center of mass and thus locomote (others include height variations and angular momentum variations that separate the ZMP and centroidal moment pivot).

Pages of this website are under the CC-BY 4.0 license.