text:
Scientific Potential of DECIGO Pathfinder and Testing GR with Space-Borne Gravitational Wave Interferometers

# Scientific Potential of DECIGO Pathfinder and Testing GR with Space-Borne Gravitational Wave Interferometers

## Abstract

DECIGO Pathfinder (DPF) has an ability to detect gravitational waves from galactic intermediate-mass black hole binaries. If the signal is detected, it would be possible to determine parameters of the binary components. Furthermore, by using future space-borne gravitational wave interferometers, it would be possible to test alternative theories of gravity in the strong field regime. In this review article, we first explain how the detectors like DPF and DECIGO/BBO work and discuss the expected event rates. Then, we review how the observed gravitational waveforms from precessing compact binaries with slightly eccentric orbits can be calculated both in general relativity and in alternative theories of gravity. For the latter, we focus on Brans-Dicke and massive gravity theories. After reviewing these theories, we show the results of the parameter estimation with DPF using the Fisher analysis. We also discuss a possible joint search of DPF and ground-based interferometers. Then, we show the results of testing alternative theories of gravity using future space-borne interferometers. DECIGO/BBO would be able to place 4–5 orders of magnitude stronger constraint on Brans-Dicke theory than the solar system experiment. This is still 1–2 orders of magnitude stronger than the future solar system mission such as ASTROD I. On the other hand, LISA should be able to put 4 orders of magnitude more stringent constraint on the mass of the graviton than the current solar system bound. DPF may be able to place comparable constraint on the massive gravity theories as the solar system bound. We also discuss the prospects of using eLISA and ASTROD-GW in testing alternative theories of gravity. The bounds using eLISA are similar to the LISA ones, but ASTROD-GW performs the best in constraining massive gravity theories among all the gravitational wave detectors considered in this article.

Gravitational Waves, DECIGO, DECIGO Pathfinder, Brans-Dicke, Massive Gravity
\catchline{history}\ccode

PACS numbers:

## 1 Introduction

### 1.1 Gravitational Waves in General

In general relativity (GR), gravitational fluctuations propagate as “ripples” in the spacetime at the speed of light, known as gravitational waves (GWs). Since gravity is weak compared to other fundamental forces, the interaction of GWs with other matters is very faint. This means that GWs can escape from highly dense or optically thick regions that cannot be probed by electromagnetic (EM) waves (e.g. very early universe before the last scattering surface). Therefore, GWs have potentials to open up a novel astronomy and cosmology.

The existence of GWs has been proved indirectly by measuring the orbital decay rate of binary pulsars, caused by the energy loss through gravitational radiation. This can be calculated analytically within the framework of GR, and this matches with the measurements beautifully [1, 2]. However, this indirect GW measurement only detects the radiated energy via GWs and not the perturbation of the spacetime itself. Therefore, the direct detection of GWs has been awaited for a long time to detect this latter effect. Currently, several ground-based detectors that are aiming for the direct detection are shifting from initial phases (first generation) to advanced phases (second-generation) such as adv. LIGO [3, 4], adv. VIRGO [5, 6], GEO-HF [7, 8] and KAGRA (formerly called LCGT) [9, 10]. These are GW interferometers that consist of two arms. When GWs pass through this detector, the length of each arm changes, modifying the interferometric pattern. KAGRA is novel and sometimes called 2.5th-generation detector since it will be cryogenically-cooled and is buried underground to reduce the thermal and seismic noises. These detectors have their best sensitivities at around 100–1000Hz. Basically, the high-frequency parts of these detectors are limited by the shot noises while the low-frequency parts are limited by the radiation pressure and the seismic noises.

Sources and sciences from GW interferometers have been summarized in Refs. [11, 12]. Promising sources for these second-generation detectors are signals from compact binaries (see Secs. 3.4 and 3.5 of Ref. [12] for a theoretical overview and Refs. [13, 14, 15, 16] for observational results) whose event rates are estimated as /year [17]. (See Ref. [18] for the review on the electro-magnetic counterparts of GW sources.) Binary signals are divided into three phases, inspiral, merger and ringdown. During the inspiral phase, the binary separation gradually shrinks due to the energy loss via GW radiation. The inspiral waveform has been studied well under the post-Newtonian (PN) formalism [19, 20, 21]. Merger phase is highly non-linear but there has been a great progress thanks to numerical relativity [22, 23, 24]. Ringdown phases are also well-studied using black hole (BH) perturbation method [25]. GWs from a merging compact binary is often called chirp signals since their frequencies get higher and higher as they head towards merger. Since the inspiral waveform modeling has been performed well using PN approach [19], it can be used to estimate binary parameters [26] and to probe cosmology as standard sirens [27]. Also, neutron-star (NS) binaries can be used to constrain equation of state (EoS) of NSs [28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. Other interesting sources are burst signals (see Sec. 3.2 of Ref. [12] for a theoretical review and Refs. [38, 39, 40, 41] for the current observational results) including supernovae [42], magnetars and cosmic strings [43, 44], and continuous signals (see Sec. 3.3 of Ref. [12] for a theoretical overview and Refs. [45, 46, 47] for observational results) from e.g. (newly-born) neutron stars. There are also stochastic GW backgrounds (see Sec. 3.6 of Ref. [12] for a theoretical review and Ref. [48] for observational results) as potential candidates for GW sources such as the ones from cosmic-strings [49], pre-big-bang models [50] and preheatings [51]. GW backgrounds associated with inflation may be detected by ground-based detectors if the spectrum has a sufficient blue tilt. There is also a proposed third-generation project called Einstein Telescope (ET) [52, 53]. The sensitivity is increased roughly by one order of magnitude compared to the second-generation ones and ET has interesting sciences which have been summarized in Ref. [54].

### 1.2 Space-Borne Interferometers

Ground-based interferometers have difficulty in detecting GW signals lower than 1–10Hz due to the seismic noise. To overcome this problem, space-borne interferometers have been proposed. Among them, the Laser Interferometer Space Antenna (LISA) has been proposed by ESA and NASA with its optimal sensitivity around GW frequency of 1mHz [55, 56]. Classic LISA consists of three spacecrafts from which laser is emitted to one another. This triangular-shaped interferometer has armlengths of km. Due to this extremely long armlength, LISA is a transponder-type (rather than a Michelson or a Fabry-Perot-type) interferometer. (I.e. Instead of reflecting the incoming laser directly, each spacecraft once receives it and emits the laser with corresponding phase.) The triangular cluster follows the same orbit as the Earth, keeping its position 20 behind the latter. The expected targets are supermassive black hole (SMBH) binaries and white dwarf (WD) binaries. The former detection may help in clarifying the mechanism of how SMBHs are formed. A system consisting of a small compact object orbiting a SMBH is called an extreme mass ratio inspiral (EMRI) [57], and GWs from this system encodes the information of the structure of the strong-field spacetime [58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76]. Unfortunately, GWs from WD binaries may mask [77, 78] other signals including primordial GW background [79, 80]. As a prototype mission, LISA path finder (LPF) will be launched before LISA [81].

Although, the classic LISA was planned and developed by NASA and ESA, the collaboration between these 2 agencies terminated in 2011. The new version of LISA (called eLISA) led by ESA alone is now under consideration. The arm-length would be shortened by 5 times compared to the classic one and the number of the arms would be reduced to 2. The satellites are slowly drifting away from Earth to save propellant. The expected operation period is 2-5 yrs. See Refs. [82, 83, 84, 85, 86, 87] for expected sources and sciences with eLISA.

Another space-borne interferometer, the Deci-Hertz Interferometer Gravitational Wave Observatory (DECIGO) has been proposed by Seto, Kawamura and Nakamura [89]. It consists of four clusters, each having a triangular interferometer. The major differences between DECIGO and LISA are that the armlengths are shortened to km and the former is a Fabry-Perot-type interferometer whereas the latter is a transponder-type as previously mentioned. Two out of four clusters are situated on the same site so that the correlation analysis can be performed to detect primordial GW background, while the other two are placed far apart so as to increase the angular resolutions of the sources [90, 91, 92] (see Fig. 1). It is most sensitive at around 0.1–1Hz. Similar interferometer, the Big Bang Observatory (BBO) has been suggested as a follow-on mission to LISA [93, 94]. Since it is a transponder-type interferometer, it can be considered as a mini-LISA with its configuration similar to the one of DECIGO. (The noise curves of adv. LIGO, ET, LISA and DECIGO/BBO are shown in Fig. 2 1.)

Since the WD/WD binary signals has a high cutoff frequency at Hz [78], the main target of DECIGO/BBO is the primordial GW background (PGWB) [89, 95, 96]. NS/NS foreground signals might mask PGWB instead of WD/WD signals, but it is likely that BBO would have enough sensitivity to subtract sufficient amount [97, 98]. See Refs. [89, 94, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114] for other cosmological and astrophysical scientific potentials of DECIGO and BBO.

As for the prototype mission, DECIGO Pathfinder (DPF) [115, 91, 116] is hoped to be launched in 2016–2017. The main goals are to test the key technologies for the space mission and to carry out observations of GWs and Earth gravity. The first space-borne GW detector, SWIM, has already been launched in 2009  [117, 116]. It is a torsion-bar type space antenna. The ground-based torsion-bar antennae have been used to place upper limits on the energy density of stochastic GWs at 0.1–1Hz  [118, 119]. The improved version of the torsion-bar antenna (TOBA) has been proposed by Ando et al[117].

DPF has a sensitivity of at . The main GW sources for DPF are the intermediate mass BH (IMBH) binaries. There are evidences that stellar-mass BHs and SMBHs exist, but there is no unambiguous detection of individual IMBHs. IMBH detection may reveal the formation mechanism of SMBHs. DPF has enough sensitivity to detect GWs from galactic IMBH binaries that might exist at the centers of globular clusters and massive young clusters (see Refs. [120, 121] for reviews on IMBHs and Ref. [122] for the possibility of IMBH binary formation). For the observation of the Earth gravity, DPF has an ability to perform complementary operation compared to other missions such as CHAMP [123], GRACE [124] and GOCE [125].

ASTROD-GW (Astrodynamical Space Test of Relativity using Optical Devices optimized for GW detection) is another space-based GW interferometer mission [126]. It consists of 3 satellites that are placed at near Lagrange points L3, L4 and L5, respectively, forming a large triangluar interferometer similar to LISA. Since its armlength is 52 times larger than that of LISA, it is more sensitive than LISA at lower frequency.

### 1.3 Testing Alternative Theories of Gravity using GWs

Among all the sciences that have been mentioned above, there is yet another very interesting and fundamentally important science that can be best probed by using GWs. It is the tests on alternative theories of gravity [127, 128, 129]. In order to solve problems like dark energy [130], dark matter and inflation [131, 132] within the context of GR, usually we need to introduce unknown matters or fields, but there are possibilities that these problems can be explained naturally by modifying gravitational theories. Also, if the classical gravitational theory is to be realized at the low energy limit of more fundamental theory like superstring theory [133, 134], the classical one does not necessarily reduce to GR due to the additional fields (e.g. dilatons) that couple to gravity.

The tests on alternative theories of gravity have been performed with great accuracies under the solar system experiments and binary pulsar observations [128, 2] and no deviation from GR has been reported so far. The former can probe only the weak field limit of the theory, while the latter can probe the strong field gravity [135] such as the effacement principle and the scalarization [136, 137, 138] in scalar-tensor theory. However, since the typical velocity of the binary component is , the system is not so dynamical. We are here interested in testing GR in the strong and dynamical field regime. The best way to perform these tests in the near future is to use GWs (especially GWs from compact binaries) since they directly contain gravitational information in the strong-field regions. Also, GWs lead to remarkably high precision GR tests in the following sense. Usually, GWs are so weak that they are buried under the detector noises. In order to dig them out from the noises, we need to perform the matched filtering analysis [139, 26]. This is nothing but taking correlations between the observed signals and the template waveforms so that if the former match with the latter, we can extract the signals out of the noises. Therefore, parameter estimations using this matched filtering technique are sensitive to deviations in GW phase since the correlations are remarkably reduced even if the phases between GWs and the templates are different only slightly. Now, it is often the case that alternative theories of gravity modify the phases of GWs from compact binaries (see Fig. 3). Let us say that we observe a signal at a frequency Hz for observation period yr. Then, the number of GW cycles (the number of phases) that we see is roughly . This shows that we would be able to detect the deviations from GR if they modify the phase at least by ! There are many works calculating how accurately we can probe alternative theories of gravity using GWs. (See e.g. Ref. [140] and references therein, and see Refs. [141, 142, 143, 144, 140, 145, 146, 147, 148, 149, 150, 151, 152] for model-independent tests. See Refs. [153, 154, 155, 156, 157] for possible tests on GR using ringdown signals.)

Two of the most important characteristics of GR are (i) it has only two tensor gravitational degrees of freedom and (ii) the graviton is massless (in other words, GWs propagate at the speed of light). In this review, we focus on simple extensions of each of the above points by (i) introducing an additional scalar field and (ii) giving a finite mass to graviton.

As for the first point, we consider Brans-Dicke (BD) theory [158] which is the simplest representative of the scalar-tensor theory [159]. This theory is parameterized by the so-called Brans-Dicke parameter which gives the (inverse) coupling between the additional (massless) scalar field and the matter fields. If we take the decoupling limit of , this theory reduces to GR. The current strongest constraint has been obtained from the solar system experiment by Shapiro time delay measurements using the Cassini satellite [160] as . There is a future space mission called ASTROD which is expected to be launched in 2021 [161]. It would perform very precise tests of GR in the solar system and the constraint on would be 3 orders of magnitude stronger than the Cassini bound. One of the remarkable points about considering gravitational radiation in this theory is that due to the additional scalar field, there exists scalar dipole radiation [162, 163, 164]! Because of this additional energy release, the binary evolution is modified compared to GR, which further makes a change in GW phase. The orbital decay of the binary pulsar PSR J1738-0333 has put a slightly weaker constraint than the solar system one [165]. There have been several works calculating the possible bounds on using future GW interferometers [166, 167, 168, 169, 170]. (See Refs. [171, 172] for the constraint on massive BD theory, Ref. [173] for the one on generic scalar-tensor theories and Ref. [174] for the bound on generic dipolar gravitational radiation.) Unfortunately, these results show that the bounds from future GW observations using adv.LIGO and LISA would not be able to exceed the one from the solar system experiment. In this review, we will show that DECIGO/BBO would be able to place 4 orders of magnitude stronger constraint than the Cassini bound [175].

For the second point, there exist many theories in which graviton acquires a finite mass  [176, 177, 178, 179, 180, 181, 182]. (See e.g. Refs. [183, 184] for recent reviews.) These theories are called massive gravity or massive graviton (MG) theories. From a pioneering work by Fierz and Pauli [176] in 1939, most of massive gravity theories have suffered from problems like Boulware-Deser ghost modes appearing in the curved background [185] or violating the Lorentz invariance [178, 179]. Only recently a self-consistent massive gravity is proposed [181, 182]. Independent of these specific massive gravity theories, the gravitational potential is modified to the one of Yukawa type. This means that the effective gravitational constant depends on the distance from the source, which modifies the Kepler’s third law. This has been tested in the solar system and the constraint has been put on the graviton Compton wavelength as cm  [186]. As for GWs, the propagation speed is modified from the speed of light which again makes a change in GW phase. GWs can also perform a model-independent tests on the graviton mass. Compared to the solar system experiment, ground-based detectors can only put comparable (or slightly stronger) constraints [187, 188, 189] whereas LISA bounds should be stronger by four orders of magnitude [187, 168, 169, 188, 190, 170, 189, 191, 145, 192, 147, 193]! This is because the deviation in the propagation speed of GWs is larger for lower frequency GWs. DECIGO/BBO would be able to place slightly weaker constraint than LISA [175] and DPF may be able to place comparable constraint to the solar system bound [194].

### 1.4 Organization

In this review, we will first explain how accurately DPF would be able to determine binary parameters by detecting GWs from galactic IMBH binaries. Then, we will describe the ability of testing GR using space-borne GW interferometers, especially DECIGO/BBO and DPF. This article is organized as follows. In Secs. 2 and 3, we explain the basic concepts, noise curves and event rates of DPF and DECIGO/BBO, respectively. In Sec. 4, we derive GWs from compact binaries, starting from the quadrupole, leading contribution for the binary with a circular orbit, and extending it to include higher PN terms and a slightly eccentric orbit. Then, in Sec. 5, we explain how to construct the observed gravitational waveforms for the spin-aligned and precessing binaries. We also introduce the inspiral-merger-ringdown hybrid waveform. In Secs. 6 and 7, we consider BD and MG theories, respectively. We explain how the GW phase is modified from GR, and review current constraints. In Sec. 8, we show the results of the binary parameter estimation using DPF, and in Sec. 9, we show the proposed constraints on BD and MG theories (and also on other theories) using future GW interferometers. Finally, we conclude in Sec. 10.

## 2 DECIGO Pathfinder

### 2.1 Basic Designs

DPF [115, 91] is a prototype mission of DECIGO to test the advanced technologies of GW space mission such as (i) a precise position measuring system with Fabry-Perot (FP) cavity, (ii) a highly stabilized laser source and (iii) the drag-free control system which shields external forces due to solar radiation and residual atmosphere. The weight of the satellite is about 350kg and it will be orbiting the earth at the Keplerian velocity with an altitude of 500km. The expected observation period is 1 year. It consists of 1-arm interferometer with armlength 30cm and the laser power of 100mW. FP cavity has not been tested in space up to now and DPF is expected to have better sensitivity than LPF which uses a Mach-Zender interferometer [81]. DPF provides new possibilities for a precise measurement of position and high-stabilized laser in space.

Each mirror is placed inside a module called housing. The relative positions of these mirrors and the frame will be measured by the local sensor and will be fed back to the satellite position using thrusters. LPF will operate at the Lagrange 1 (L1) point where the gravitational environment is stable, while DPF will demonstrate it in an earth orbit. This will open a new window for future space missions.

DPF GW observations aim at the frequency of 0.1–1Hz which is important because the ground-based interferometers and LPF are not sensitive in this frequency range. Also DPF data are expected to be more complicated than the one from ground-based interferometers due to the satellite orbital motion and the effects of the earth, and hence the development of data analysis technique for DPF has a significant meaning.

### 2.2 Noise Spectrum and Observable Range of DPF

The detected signal can be expressed as the combination of the GW signal and the noise . If the noise is stationary, the noise spectral density can be defined as

 ⟨~n∗(f)~n(f)⟩=δ(f−f′)12Sn(f), (1)

where represents the Fourier component of the noise and the angle brackets denote the expectation value. The root noise spectral density of DPF is given as

 Sn(f)=1.0×10−30(2π)4(f1Hz)−4+4.0×10−30 Hz−1, (2)

where the first term corresponds to the acceleration noises while the second term represents the laser frequency noise. It is shown in Fig. 4 as the (red) thick solid curve. We set the lower frequency cutoff at Hz due to the Earth gravity. It may be possible to improve the sensitivity in which the higher frequency part is now limited by the shot noise. We call this the “adv. DPF” whose noise spectral density is given by [195]

It is shown in Fig. 4 as the (green) thin solid curve. We also show the lower cutoff frequencies of the adv. LIGO (magenta thick dotted) and KAGRA (blue thin dotted) 2.

In Fig. 5, we show the (sky-averaged) observable range of DPF and adv. DPF as the (red) thick solid curve and the (green) thin solid curve, respectively. We also show the ones of adv. LIGO and KAGRA with the same curves as in Fig. 4. The (black) dashed horizontal line at kpc corresponds to the galactic center. We also plot the possible GW sources at the centers of the globular clusters shown in Table 2.2, assuming that they consist of equal-mass IMBH binaries.

### 2.3 Beam Pattern Functions and the Effect of Detector Motion

The beam pattern functions of 1-armed interferometer are given by [194]

 F+(θ,ψ) = 12sin2θcos2ψ, (4) F×(θ,ψ) = 12sin2θsin2ψ. (5)

Here, is the angle between the arm and the incoming GW, and represents the polarization angle (see Fig. 6). The sky-averaged values of them are . In Fig. 4, we also show the sky-averaged GW amplitudes of a equal-mass IMBH binary in Centauri (light blue thick dotted-dashed) and a equal-mass one in NGC 6752 (black thin dotted-dashed). We assume the dimensionless spin parameter of these BHs to be whose value is predicted to be the most probable [206] (see also Ref. [207]). We use the phenomenological inspiral-merger-ringdown hybrid waveforms [208, 209] (see Sec. 5.3) to estimate these amplitudes. One can see that GW frequency of IMBH binary in Centauri is too low for the ground-based interferometers, and hence its GW signals become unique sources for DPF. On the other hand, the one in NGC 6752 can be detected with both DPF and the ground-based ones. This suggests that it may be possible to perform joint searches between these detectors [194], which we explain in Sec. 8.1.

For later use, we introduce the celestial coordinates shown in Fig. 7. -axis points the vernal equinox, while -axis points the north celestial pole and is orthogonal to the celestial plane. The source direction and the polarization angle can be re-expressed in terms of the source direction and the direction of the orbital angular momentum as [194]

 cosθ = cosφD(t)siniDcos¯θS (6) −[cos{φE(t)−¯ϕS}cosφD(t)cosiD+sin{φE(t)−¯ϕS}sinφD(t)]sin¯θS. cosψ = [cos¯θS{−cos{φE(t)−¯ϕL}sinφD(t)+sin{φE(t)−¯ϕL}cosφD(t)cosiD}sin¯θL (7) +{cos{φE(t)−¯ϕS}sinφD(t)cos¯θL −cosφD(t){sin{φE(t)−¯ϕS}cos¯θLcosiD+sin{¯ϕL−¯ϕS}sin¯θLsiniD}}sin¯θS] [(1−cos2θ){1−(sin¯θSsin¯θLcos(¯ϕL−¯ϕS)+cos¯θScos¯θL)2}]−1/2.

### 2.4 Event Rate

In this subsection, we estimate the event rate for (I) Equal-mass BH binaries and (II) IMRIs using DPF.

#### Equal-mass BH Binaries

DPF can detect GW signals from galactic IMBH binaries. Two of the candidates where IMBH binaries might exist in our Galaxy are at the centers of globular clusters (GCs) and galactic massive young clusters (GMYCs).

For the former case, we select 21 globular clusters out of 150 globular clusters that have been found. We determine the total mass of the possible IMBH binary by applying the formula obtained by Tremaine et al. [210] as

 M=1.35×108(σ200km/s)4M⊙. (8)

We assume the mass ratio of IMBH binary to be 1. The results are listed in Table 2 of Ref. [194]. We plot IMBH binary of each globular cluster in Fig. 5 and count how many of them lie within the observable range of DPF and Adv. DPF. One sees that 2 out of 26 globular clusters (5 shown in Table 2.2 + 21 shown in Table 2 of Ref. [194]) might contain IMBH binaries detectable with DPF on average. Then, the detection rate of the IMBH binaries in globular clusters is given by

 (9)

Here, following Ref. [211], we divide the number of globular clusters by the age of the universe. This is because only one IMBH binary can be formed over its lifetime in each cluster. When we use the adv. DPF, 13 out of 26 IMBH binaries in the galactic globular clusters lie withing the observable reach of the detector. This makes the event rate larger than Eq. (9) by a factor .

For the latter case, currently, more than 10 galactic massive young clusters (GMYCs) have been discovered [212]. Grkan et al. [122] performed numerical simulations and found that two IMBHs may form in GMYCs if the initial binary fraction is relatively large. After IMBH binary formation, it shrinks due to the dynamical friction with the cluster stars. The timescale of the dynamical encounters would be Gyr  [211]. For simplicity, we assume that IMBH binaries are all located at the galactic center. From Fig. 5, one sees that DPF is not sensitive enough to detect GW signals from IMBH binaries in GMYCs. On the other hand, when we use adv. DPF, it has ability to detect all of the IMBH binaries in GMYCs. Following Fregeau et al. [211], we assume that the number of star clusters massive enough to form IMBH binaries is the same as the one of globular clusters, and 10 of them actually produce IMBH binaries. In our galaxy, about 150 globular clusters [201] have been observed. This means that at least 15 IMBH binaries are expected to be within the reach of adv. DPF. The detection rate of the IMBH binaries in GMYCs with adv. DPF can be roughly estimated as .

#### IMRIs

After supernova explosions, many stellar-mass BHs of would form. These would sink to the cluster cores due to the mass segregation, and their number density is expected to be comparable to the one of the main sequence stars, typically pc. In this subsection, we consider IMBH binaries with intermediate-mass ratio inspirals (IMRIs). In Fig. 8, we show observable ranges of DPF (red thick solid) and adv. DPF (green thin solid) with . One sees that DPF is not sensitive enough to detect any GWs from IMRIs with , but adv. DPF has ability to detect these signals from some of GCs and the galactic center.

For GMYCs, from Fig. 8, one sees that even adv. DPF is not sensitive, on average, to an IMRI signal of at the galactic center. We found that 1.7 of these binaries have SNRs greater than 5 due to the optimal orientations [194]. On the other hand, the merger rate due to three-body interaction per cluster becomes  [213] 3. We make a conservative assumption that 50 GMYCs may contain IMBHs (see e.g. Refs. [212, 213]). In this case, the total detection rate becomes

 ˙NGMYC,tot=0.017×50×2×10−8yr−1=1.7×10−8yr−1. (10)

For GCs, we found that the merger rate is limited by the requirement that the merger rate in a single cluster cannot exceed  [206, 194], where is the current age of the cluster. Otherwise, the mass of the larger BH has been increased considerably. From Fig. 8, one sees that IMRI signals from about 3 out of 26 globular clusters might be detected by adv. DPF. This leads to the estimate that it can detect IMRI signals from about globular clusters in total, and hence the detection rate of globular clusters in total becomes

 ˙NGC,tot≲1.2×10−7yr−1. (11)

## 3 Decigo/bbo

### 3.1 Basic Designs

DECIGO is first proposed by Seto, Kawamura and Nakamura [89]. It consists of four triangular sets of detectors whose configuration is shown in Fig. 1. This effectively corresponds to eight individual interferometers. Each triangular detector has an arm-length of km. Its primary goal is to detect PGWB with the GW energy density of . Since the WD/WD confusion noise will have a cutoff frequency at around 0.2 Hz  [78], DECIGO has an advantage on detecting this source over LISA (see Refs. [97, 98, 214] for the discussions of the NS/NS confusion noise). Therefore two of the four triangular detectors are located on the same site forming a star of David so that correlation analysis [215, 216] can be performed to detect PGWB. The rest of the two detectors are placed far apart to increase the angular resolutions of the source locations.

BBO has almost the same constellation as DECIGO. The main difference is that while DECIGO is a Fabry-Perot type interferometer, BBO is a transponder-type interferometer. BBO has arm-lengths of km.

### 3.2 Noise Spectrum

The noise spectrum of BBO is given as follows. The non-sky-averaged instrumental noise spectral density for BBO is obtained from Ref. [94] as [88]

 Sinstn(f)=[1.8×10−49(f1Hz)2+2.9×10−49+9.2×10−52(f1Hz)−4 ]Hz−1. (12)

It has 20/3 times better sensitivity than the one for the sky-averaged sensitivity [169]. Apart from instrumental noise, there are astrophysical foreground confusion noises. These confusion noise spectral densities and the energy densities of gravitational waves are related as [97]

 Sconfn=35πf−3ρcΩGW. (13)

Here, is the critical energy density of the Universe and

 ΩGW≡1ρcdρGWdlnf (14)

is the energy density of GWs per log frequency normalized by . The energy density of GWs that originate from extra-galactic WD binaries has been estimated as which leads to the noise spectral density of [78]

 Sex−galn(f)=4.2×10−47(f1 Hz)−7/3 Hz−1. (15)

On the other hand, the one from galactic WD binaries has been calculated as [77]

 (16)

We multiply and by a factor , which corresponds to the high frequency cutoff for the white dwarf confusion noises. We also have to take into account the confusion noise from NS binaries. Its noise spectral density is estimated as [97, 214]

 SNSn(f)=1.3×10−48(f1 Hz)−7/3(˙n010−7 Mpc−3yr−1)Hz−1, (17)

where denotes current merger rate density of NS/NS binaries. Putting altogether, the total noise spectral density for BBO becomes

 Sn(f) = (18) +Sex−galn(f)F(f)+FcleanSNSn(f).

Here is the number density of galactic white dwarf binaries per unit frequency, which is given by [217]

 dNdf=2×10−3Hz−1(f1 Hz)−11/3. (19)

is the average number of frequency bins that are lost when each galactic binary is fitted out. The factor in front of represents our assumption of the fraction of GWs that cannot be removed after foreground subtraction. In this review, we assume that NS/NS foregrounds can be subtracted down to the level below the instrumental sensitivity. The lower and higher frequency ends of the BBO sensitivity band are set as Hz and Hz, respectively. The noise spectrum of BBO is shown as a red thick solid curve in Fig. 2.

DECIGO has been proposed with 3-4 times less sensitive spectrum than BBO. Its instrumental noise spectrum is given by

 Sinstn,DECIGO(f)=5.3×10−48[(1+y2)+2.3×10−7y4(1+y2)+2.6×10−8y4]Hz−1, (20)

where with Hz. However, this is not the fixed design sensitivity and there is a project going on to improve the sensitivity to the same level as BBO. The frequency range of DECIGO/BBO is given as

 (flow,fhigh)=(10−3Hz,102Hz). (21)

### 3.3 Event Rate

The promising sources for DECIGO/BBO are inspirals of NS binaries. In this subsection, following Cutler and Harms [97], we derive the detection rate of these sources. Here, we adopt the model where NS/NS binaries only exist at redshift below . (This is inferred from observations. See the discussion below.) Since DECIGO/BBO has enough sensitivity to detect the farthest NS/NS binaries considered here, the merger rate corresponds to the detection rate. The NS/NS merger rate can be written as

 ˙NNSNS=∫∞04π[a0r1(z)]2˙n(z)dτ1dzdz, (22)

where

 a0r(z)=∫z0dz′H(z),dτdz=1(1+z)H(z). (23)

Here, represents the Hubble parameter at redshift . For -CDM cosmology, it is given by

 H(z)≡H0√(1+z)3Ωm+ΩΛ. (24)

is the merger rate of NS/NS binaries at redshift and it can be re-expressed in terms of the current merger rate and the redshift dependence as . For , we adopt the following piecewise linear fit  [97, 218] based on observations discussed in Ref. [219];

 R(z)=⎧⎪ ⎪⎨⎪ ⎪⎩1+2z(z≤1)34(5−z)(1≤z≤5)0(z≥5). (25)

Then, we finally obtain the detection rate as

 ˙NNSNS=105(˙n010−7Mpc−3yr−1)yr−1. (26)

The detection rate for BH/NS binaries are expected to be roughly 1/10 of NS/NS detection rate.

We can also estimate the NS/NS foreground as follows. First, we apply the convenient formula given by Phinney [220] as

 ΩNSGW = 8π5/391H20M5/3f2/3∫∞0dz˙n(z)(1+z)4/3H(z) (27) = 3.74×10−12h−372(M1.22M⊙)5/3(f1Hz)2/3(˙n010−7Mpc−3yr−1).

By using Eqs. (13) and (27), we obtain

 √SNSn=1.76×10−24h72(M1.22M⊙)5/6(˙n010−7Mpc−3yr−1)1/2(f1Hz)−7/6. (28)

## 4 Gravitational Waves from Compact Binaries

In this section, we focus on the GWs from inspiral binaries composed of compact stars such as BHs and NSs within the context of GR. In Sec. 4.1, we first derive the GWs from a circular binary at the leading order. Then, In Sec. 4.2, we extend the result to higher order terms. Finally, in Sec. 4.3, we consider binaries with small eccentricity.

In this section, we derive the leading order gravitational waves from compact binaries, namely the quadrupole radiation. The binary circularizes quickly during its orbital decay due to the gravitational radiation reaction, and hence we restrict our attention to the binary with circular orbit in this section. Also, we assume that the velocity of each body is sufficiently small compared to the speed of light.

Let us assume that we have two compact objects with masses and in the - plane with orbital separation at a distance from the detector. We define the total and the reduced mass as and , respectively. The binary components obey the Kepler’s Law at the leading order with orbital angular velocity defined as

 Ω≡√Mta3. (29)

As shown in Fig. 9, we denote the positions of these objects in the center of mass frame as and , which can be expressed as

 (x1,y1) = (m2Mtacos(Ωt+π2),m2Mtasin(Ωt+π2)), (30) (x2,y2) = (31)

By taking the time derivative of the quadrupole moment twice, we obtain

 h+ = A+cos2Ωt, (32) h× = A×sin2Ωt. (33)

Here, we have defined

where is the unit orbital angular momentum vector and is the unit vector pointing towards the direction of GW propagation. From this result, we see that the monochromatic GW frequency can be written as at the leading order.

Next, we take the effect of radiation reaction into account and estimate the GW phase  [26]. First, the total energy of the binary system is given as

 E=−μMt2a. (36)

Next, by using the quadrupole formula

 dEGWdt=15⟨...Mij...Mij−13(...Mkk)2⟩, (37)

we have the radiated energy due to GWs as

 dEGWdt=325μ2M3ta5. (38)

By using the balancing equation

 dEdt=−dEGWdt, (39)

the evolution of the binary separation can be expressed as

 ˙a=−aEdEdt=−645μM2ta3, (40)

which can be solved to yield

 a(t)=(2565μM2t(t0−t))1/4, (41)

where is the coalescence time. This can be turned into the frequency evolution as

 ˙f=˙Ωπ = −32M1/2tπ˙aa5/2 (42) = 965π8/3M5/3f11/3,

where is the chirp mass. This equation can be solved as

 f=(5256)3/81πM5/81(t0−t)3/8. (43)

By integrating this by once, we get the phase as

 ϕ(t)=∫2πfdt=−2(15M−1(t0−t))5/8+ϕ0, (44)

with representing the coalescence phase.

Next, we derive the waveform in the Fourier domain defined as

 ~h(f)≡∫∞−∞e2πifth(t)dt. (45)

Since we have assumed that the velocity of each binary component is much smaller than the speed of light, the amplitude and the phase satisfy the following conditions: and . Under this situation, we can apply the stationary phase approximation, where the Fourier component of a function becomes [26]

 ~B(f)≈12A(t)(dfdt)−1/2expi(2πft−ϕ(f)−π/4). (46)

Here, is the time that satisfies and . From Eqs. (43) and (44), we have

 t(f) = t0−5M(8πMf)−8/3, (47) ϕ(f) = ϕ0−2(8πMf)−5/3. (48)

Therefore, we obtain

 ~h+(f) = Af−7/6(1+(^L⋅n)22)eiΨ+(f), (49) ~h×(f) = Af−7/6(^L⋅n)eiΨ×(f), (50)

where the amplitude and phases are given as

 A = √5241π2/3M5/6DL, (51) Ψ+(f) = 2πft0−ϕ0−π4+Ψ0PN(f), (52) Ψ×(f) = Ψ+(f)+π2, (53)

respectively. Here, is defined as

 Ψ0PN(f)≡3128(πMf)−5/3, (54)

and is the luminosity distance given as

 DL=(1+z)∫z0dz′H(z′). (55)

### 4.2 Higher Post-Newtonian Corrections

In the previous section, we derived the waveform at the leading order, where we used the quadrupole formula for the radiated energy flux and assumed that the binary orbit is of Newtonian. In this section, we extend the previous analysis to higher orders. Since we go beyond Newtonian, this expansion is called the post-Newtonian (PN). Our expansion parameter is where is the typical velocity of the binary component, and the terms that are proportional to are of PN orders. In this section, we consider the extension up to 2PN order [221].

Furthermore, we neglect the PN corrections to the amplitude. This is because when we take the correlation between two different waveforms, it is more sensitive to the deviation in the phase than the one in the amplitude. This type of waveform is called the restricted PN waveform.

In Ref. [221], expressions for , and are shown to 2PN order. By using them, up to this PN order, becomes

 ˙f = 965π8/3M5/3f11/3[1−(743336+114η)x+(4π−β)x1/2 (56) +(3410318144+136612016η+5918η2+σ)x2],

where is the symmetric mass ratio with PN expansion parameter . The spin-orbit and the spin-spin couplings are given as

 β = 112∑i=1,2(113m2iMt+75η)^L⋅χi, (57) σ = 148η{−247χ1⋅χ2+721(^L⋅χ1)(^L⋅χ2)}, (58)

where the dimensionless spin parameter is defined as with denoting the spin angular momentum of the -th compact object. The gravitational wave phase in the Fourier domain can be calculated as

 Ψ+(f) = (59) −4(4π−β)x3/2+(15293365508032+27145504η+308572η2−10σ)x2]. Ψ×(f) = Ψ+(f)+π2. (60)

### 4.3 Gravitational Waves from Binaries with Slightly Eccentric Orbit

In this section, we consider gravitational waves from a binary that have a small eccentricity . We can express the eccentric orbit as

 rorb=a(1−e2)1+ecosΩt, (61)

where the orbital angular velocity can be written as

 Ω=√Mta(1−e2)r2. (62)

By using the quadrupole formula, the radiated energy flux can be estimated as

 dEGWdt=325η2M5ta5(1−e2)7/2(1+7324e2+3796e4), (63)

and similarly, the radiated angular momentum flux is calculated as

 dJGWzdt = 25εijk⟨...Qxy(¨Qyy−¨Qxx)+¨Qxy(...Qxx−...Qyy)⟩ (64) = 3254η2M9/2ta7/2(1−e2)2(1+78e2).

Next, we derive the evolution of the semi-major axis and the eccentricity . Since the orbital energy of the binary is given in Eq. (36) and the angular momentum is given as

 Jz=μr2orbΩ=ηM2t√Mt√a(1−e2) (65)

at the Newtonian order, by using Eqs. (63) and (64), the evolutions of and can be calculated as

From these equations, we obtain

 daa=19121+7324e2+3796e4e(1−e2)(1+121304e2)de. (68)

Therefore, by integrating once, we obtain the relation between and as

 aai=1−e2i1−e2(eei)12/19⎡⎣1+121304e21+121304e2i⎤⎦, (69)

where and denote the initial semi-major axis and eccentricity, respectively.

Finally, we derive the gravitational waveform (in the Fourier domain) from the compact object with the small eccentric orbit [222, 97]. By substituting Eq. (69) into the GW frequency , we obtain

 f(e) = fiζ(ei)ζ(e), (70) ζ(e) ≡ e18/19(1+121304e2)13052299(1−e2)3/2, (71)

where is the frequency at . Especially, when , we have

 e≈ei(ffi)−19/18, (72)

so that up to , the asymptotic eccentricity invariant,

 Ie≡e2f19/9 (73)

is conserved.

The evolution the frequency is given by

 ˙f=965π8/3M5/3f11/3(1+15724e2)+O(e4), (74)

for . By substituting Eq. (72) into the above equation, we have

 ˙f=965π8/3M5/3f11/3(1+157


text:
# Home

Score : $100$ points

### Problem Statement

You are given a string $S$ of length exactly $9$ consisting of digits. One but all digits from 0 to 9 appear exactly once in $S$.

Print the only digit missing in $S$.

### Constraints

• $S$ is a string of length $9$ consisting of digits.
• All characters in $S$ are distinct.

### Input

Input is given from Standard Input in the following format:

$S$

### Output

Print the only digit missing in $S$.

### Sample Input 1

023456789

### Sample Output 1

1

The string 023456789 only lacks $1$. Thus, $1$ should be printed.

### Sample Input 2

459230781

### Sample Output 2

6

The string 459230781 only lacks $6$. Thus, $6$ should be printed.

Note that the digits in the string may not appear in increasing order.

format: plain


text:
# Multivariate Gaussian distribution

by Robert M. Johnson - Sat 02 July 2016
Tags:

In this article, we work out the mathematics of what happens when we consider Gaussian probability distributions in more than one dimension. Let's start. Recall the Gaussian distribution. It's probability density function is given by:

$N(x \, | \, \mu_x, \sigma_x^2) = \frac{1}{\sqrt{2 \pi \sigma_x^2}}e^{-\frac{1}{2}\frac{(x-\mu_x)^2}{\sigma_x^2}}$

$\int_{-\infty}^{\infty} \, N(x \, | \, \mu_x, \sigma_x^2) \, \text{dx} = 1$

The above Gaussian is univariate because it involves one variable, what if we wish to involve more than one variable? Let's see how we might construct this:

$N(x_1 \, | \, \mu_{x_1}, \sigma_{x_1}^2) = \frac{1}{\sqrt{2 \pi \sigma_{x_1}^2}}e^{-\frac{1}{2} \frac{(x_1-\mu_{x_1})^2}{\sigma_{x_1}^2}}$

$N(x_2 \, | \, \mu_{x_2}, \sigma_{x_2}^2) = \frac{1}{\sqrt{2 \pi \sigma_{x_2}^2}}e^{-\frac{1}{2} \frac{(x_2-\mu_{x_2})^2}{\sigma_{x_2}^2}}$

Here, $x_1$ and $x_2$ are different continous random variables. If we assume they are independent:

$N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) = N(x_1 \, | \, \mu_{x_1}, \sigma_{x_1}^2)N(x_2 \, | \, \mu_{x_2}, \sigma_{x_2}^2)$

Our new distribution remains normalized since:

$\int_{-\infty}^{\infty} N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) dx_1 dx_2 = \int_{-\infty}^{\infty} N(x_1 \, | \, \mu_{x_1}, \sigma_{x_1}^2) dx_1 \int_{-\infty}^{\infty}N(x_2 \, | \, \mu_{x_2}, \sigma_{x_2}^2) dx_2 = 1 \cdot 1 = 1$

Thus:

$N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) = \frac{1}{\sqrt{2 \pi \sigma_{x_1}^2}} \frac{1}{\sqrt{2 \pi \sigma_{x_2}^2}} e^{-\frac{1}{2} \frac{(x_1-\mu_{x_1})^2}{\sigma_{x_1}^2} -\frac{1}{2} \frac{(x_2-\mu_{x_2})^2}{\sigma_{x_2}^2}}$

We can clean this up a bit:

$N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) = \frac{1}{(2 \pi)^{\frac{D}{2}}} \frac{1}{(\sigma_{x_2}^2 \sigma_{x_1}^2)^{\frac{1}{2}}} e^{-\frac{1}{2} \sum_i^{D} \frac{(x_i-\mu_{x_i})^2}{\sigma_{x_i}^2}}$

Here, $D$ is the number of dimensions. In our current case, $D=2$.

## Getting Matrices Involved

We can clean this up some more. Consider the covariance matrix of this system:

$\Sigma = \begin{bmatrix} \sigma_{x_{11}}^2 & \sigma_{x_{12}}^2 \\ \sigma_{x_{21}}^2 & \sigma_{x_{22}}^2 \end{bmatrix}$

Well, by our independence assumption, we have no correlations! Thus:

$\Sigma = \begin{bmatrix} \sigma_{x_{11}}^2 & 0 \\ 0 & \sigma_{x_{22}}^2 \end{bmatrix}$

Nice! Let's calculate the determinant of this covariance matrix:

$\big| \, \Sigma \, \big| = \sigma_{x_{11}}^2 \sigma_{x_{22}}^2 - 0 \cdot 0 = \sigma_{x_{11}}^2 \sigma_{x_{22}}^2$

Using this, our equation now looks like:

$N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) = \frac{1}{(2 \pi)^{\frac{D}{2}}} \frac{1}{ | \, \Sigma \, | ^{\frac{1}{2}}} e^{-\frac{1}{2} \sum_i^{D} \frac{(x_i-\mu_{x_i})^2}{\sigma_{x_i}^2}}$

Note this looks a bit confusing with the use of two $\Sigma$ symbols. The one far to the right is stil used for summation, while the one on the left represents our covariance matrix.

We can still clean up more. Consider this term:

$\sum_i^{D} \frac{(x_i-\mu_{x_i})^2}{\sigma_{x_i}^2} = \sum_i^{D} (x_i-\mu_{x_i}) \frac{1}{\sigma_{x_i}^2} (x_i-\mu_{x_i})$

We can introduce some column vectors:

$\vec{\mu} = \begin{bmatrix} \mu_{x_1} \\ \mu_{x_2} \end{bmatrix}$ $\vec{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$

Also, note that the inverse of the covariance matrix is

$\Sigma \Sigma^{-1} = I = \begin{bmatrix} \sigma_{x_{11}}^2 & 0 \\ 0 & \sigma_{x_{22}}^2 \end{bmatrix} \begin{bmatrix} \frac{1}{\sigma_{x_{11}}^2} & 0 \\ 0 & \frac{1}{\sigma_{x_{22}}^2} \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$

Now, we just need to put all of this stuff together. Consider:

$\begin{bmatrix} x_1 - \mu_{x_1} ,\; x_2 - \mu_{x_2} \end{bmatrix} \begin{bmatrix} \frac{1}{\sigma_{x_{11}}^2} & 0 \\ 0 & \frac{1}{\sigma_{x_{22}}^2} \end{bmatrix} \begin{bmatrix} x_1 - \mu_{x_1} \\ x_2 - \mu_{x_2} \end{bmatrix}$

This is actually the same as:

$\sum_i^{D} \frac{(x_i-\mu_{x_i})^2}{\sigma_{x_i}^2}$

Thus,

$\sum_i^{D} \frac{(x_i-\mu_{x_i})^2}{\sigma_{x_i}^2} = (\vec{x} - \vec{\mu})^T \, \Sigma^{-1} \, (\vec{x} - \vec{\mu})$

## The Final Form

Finally, our equation looks like:

$N(\vec{x} \, | \, \mu, \Sigma) = N(x_1, x_2 \, | \, \mu_{x_1}, \mu_{x_2}, \sigma_{x_1}^2, \sigma_{x_2}^2) = \frac{1}{(2 \pi)^{\frac{D}{2}}} \frac{1}{ | \, \Sigma \, | ^{\frac{1}{2}}} e^{-\frac{1}{2} (\vec{x} - \vec{\mu})^T \, \Sigma^{-1} \, (\vec{x} - \vec{\mu})}$

Clearly, the following equation is much easier to write down and work with:

$N(\vec{x} \, | \, \mu, \Sigma) = \frac{1}{(2 \pi)^{\frac{D}{2}}} \frac{1}{ | \, \Sigma \, | ^{\frac{1}{2}}} e^{-\frac{1}{2} (\vec{x} - \vec{\mu})^T \, \Sigma^{-1} \, (\vec{x} - \vec{\mu})}$

Although we worked out the math for the case $D=2$, this generalizes to higher dimensions as well. Well, there you have it, the Gaussian in higher dimensions. And remember, this form assumes the $x_i$ random variables involved are all independent.


text:
# Calculate $\int_S \vec{A} \cdot \hat{n} \ dS, S=\{x^2+y^2+z^2=R^2, z \ge 0\}$

Calculate $$\int_S \vec{A} \cdot \hat{n} \ dS$$

$$\vec{A}=4 \hat{i}$$

Using spherical coordinates:

\begin{cases} \hat{r}=\sin \theta \cos \phi \ \hat{i} + \sin \theta \sin \phi \ \hat{j}+ \cos \theta \ \hat{k} \\ \hat{\theta}=\cos \theta \cos \phi \ \hat{i} + \cos \theta \sin \phi \ \hat{j}- \sin \theta \ \hat{k} \\ \hat{\phi}=-\sin \phi \ \hat{i} + \cos \phi \ \hat{j} \end{cases}

Considering:

$$\hat{n}=\hat{r}$$

$$\hat{r} \cdot \hat{i}=\sin \theta\cos \phi$$

Jacobian determinant: $$r^2 \ \sin \theta$$

The surface integral:

$$\int_S \vec{A} \cdot \hat{n} dS=\int_0^{\frac{\pi}{2}} \left( \int_0^{2 \pi} R^2 \sin \theta (4 \sin \theta \cos \phi) \ d\phi \right) d\theta$$

$$\int_0^{2 \pi} \cos \phi=0$$

Then:

$$\int_S \vec{A} \cdot \hat{n} \ dS=0$$

Is it correct?

Can I apply the Divergence theorem in this case?

Thanks!

It is correct (at least from a quick glance), and yes, you could have applied the divergence theorem by noting that $$A = \nabla f$$ where $$f(x, y, z) = 4x$$ The integral of $f$ over the solid ball (which the div. theorem would have you compute) is evidently zero, because the half-ball ($z \ge 0$) has a reflection symmetry through the plane $x = 0$, while the integrand, $f$, changes sign under this symmetry. But you'd need to take into account the integral of $A$ over the other part of the boundary of the half-ball (the disk $x^2 + y^2 = R^2, z = 0$), and see that this, too, is zero, before the job would be complete. [It's easy, because $A$ and the surface normal are orthogonal on that disk, but still...you have to do it.]


text:
# Why complex discrete Fourier transform?

1. Dec 12, 2012

### Jamal26

I've been trying to figure out why it's standard to use complex discrete Fourier transforms instead of just the real version. It's discussed a bit here.

http://dsp.stackexchange.com/questions/1406/real-discrete-fourier-transform

As far as I can tell there's a hypothetical efficiency argument. Am I wrong? I notice that a lot of computer algorithms are using the real version. For instance ffmpeg tends to use the real version as far as I can tell.

It seems to me that before computers, complex numbers were used to aid in hand computation. For instance it's easier to view sinusoidal functions as a power function of e than it is to calculate sines and cosines when you're doing it by hand.

Am I completely wrong about this analysis? In the case of Fourier transforms do I get missing data or have reduced accuracy when I use the real version instead of the complex version?

Thanks.

2. Feb 19, 2014

### AlephZero

THe stackexchange answer is pretty good. The bottom line is, there are some irritating sign-changes between $A \cos\omega t + B\sin\omega t$ and the real part of $(A + iB)e^{i\omega t}$. The math works out tidier in the general case where both the data and its transform are complex.

FWIW The stackexchange comment
is not always true. If you consder the behaviour of a rotating machine, the simple way to represent its response in the time domain is using complex variables, in a similar way that electrical engineers use "phasors".

I think the "computational efficiency argument" is entirely hypothetical. if you want efficiency (and numerical accuracy) it is not usually a good plan to make a naive translation of math formulas using complex numbers, into computer code using complex variables. But aside from that, you can use exactly the same algorithm for doing an FFT of 2n real variables or n complex variables, apart from a different interpretation of the order of the input and output data which has no effect on the computing time or cost.


text:
Thesis Open Access

# Assessment of the utility of the ESA SST CCI product for environmental research using data from meteorological buoys in the North Atlantic

Binner, Hannah

### JSON-LD (schema.org) Export

{
"inLanguage": {
"alternateName": "eng",
"@type": "Language",
"name": "English"
},
"description": "
This study assessed the accuracy and utility of the European Space Agency’s Sea Surface Temperature Climate Change Initiative (ESA SST CCI) product, which was established through the Global Ocean Data Assimilation Experiment (GODAE) and the GODAE High-Resolution SST Pilot-Project (GHRSST-PP). The product, the GHRSST Multi-product ensemble (GMPE), delivered Sea Surface Temperature (SST) measurements as a global daily median observation. These daily observations were aimed to be validated through the use of in-situ data, using four of the Met Éireann meteorological buoys. For the validation, the satellite images were processed in order to obtain the SST measurement that the satellite sensors had taken at the location of each of the four in-situ buoys and compared to the SST measurements taken by the in-situ buoys. Overall, satellite observations were found to have lower SST observations, while the buoy data showed a higher variability. Seasonal variability reflected in the satellite sensor and the buoy data. It was validated, that the satellite sensor data was accurate within 0.18 °C for the M3 buoy, 0.22 °C for the M6 buoy, 0.28 °C for the M5 buoy and 0.85 °C for the M4 buoy. It was also shown that the data was more reliable for summer months at the M3 and M5 buoy locations and least reliable for the M4 buoy location, which was attributed to a potential calibration error. Results of the M6 buoy location were attributed with a large uncertainty due to outliers.
\n\n
The utility of the SST variable GMPE data for use in Irish waters is suggested for applications in fisheries and aquaculture.
",
"creator": [
{
"affiliation": "Department of Biological, Environmental and Earth Science, University College Cork",
"@type": "Person",
"name": "Binner, Hannah"
}
],
"url": "https://zenodo.org/record/1450367",
"datePublished": "2018-03-16",
"headline": "Assessment of the utility of the ESA SST CCI product for environmental research using data from meteorological buoys in the North Atlantic",
"version": "1.0",
"keywords": [
"Sea Surface Temperature",
"Environmental Science",
"Oceanography",
"ESA CCI",
"North Atlantic",
"Satellite",
"Earth Observation",
"Fisheries and Aquaculture application"
],
"@context": "https://schema.org/",
"identifier": "https://doi.org/10.5281/zenodo.1450367",
"@id": "https://doi.org/10.5281/zenodo.1450367",
"@type": "ScholarlyArticle",
"name": "Assessment of the utility of the ESA SST CCI product for environmental research using data from meteorological buoys in the North Atlantic"
}
69
55
views


text:
# 5-digit puzzle code looking for solution

Can somebody help me solve this, or can you teach me how?

4 7 2 9 1 - One number is correct but not in right position
9 4 6 8 7 - One number is correct but not in right position
3 1 8 7 2 - Two numbers are correct but only one is in right position
1 5 7 3 9 - Two numbers are correct and both in right position

• Hi, welcome to Puzzling.SE! Where does this puzzle come from? If someone else made it, you need to provide a link back to wherever you found it. – F1Krazy Apr 12 '20 at 11:57
• its a game from my boss dont know where he get it,, so im asking someone to solve it.really got headache solving this – Ayhs Sepsep Apr 12 '20 at 12:06
• we did came up with the same answer.. hope its correct. thank u very much – Ayhs Sepsep Apr 13 '20 at 12:32
• If Oray's answer does turn out to be correct, you can mark it as correct by clicking the green tick mark next to his answer. If it isn't, but your boss tells you the correct solution, you can post that as an answer yourself. – F1Krazy Apr 13 '20 at 12:34
• still waiting for his comment on my answer – Ayhs Sepsep Apr 13 '20 at 12:49

65032

because

I believe your boss' intention was using every digit only once even though it was not given.

Otherwise, there are lots of answers.

How I solved;

Assuming

Every digit needs to be used only once

First two 5 digit numbers (guesses)

have 3 digits in common, so I eliminate them first: $$9$$,$$7$$,$$4$$. The reason is that if we choose any one these digits, we need to eliminate at least 6 numbers directly: $$1$$,$$2$$,$$6$$,$$8$$ and two more numbers among the ones I eliminated. So it should not be the case.

Then

I choose couple from the first wo digit numbers with the left digits, and at the end of trial and error I found the answer after like a minute. (such as {1,8} {1,6}, {2,8}, {2,6}).


text:
# arch.univariate.ZeroMean.starting_values¶

ZeroMean.starting_values()

Returns starting values for the mean model, often the same as the values returned from fit

Returns

sv – Starting values

Return type

ndarray


text:
# Manacher's Algorithm - Finding all sub-palindromes in $O(N)$

## Statement

Given string $s$ with length $n$. Find all the pairs $(i, j)$ such that substring $s[i\dots j]$ is a palindrome. String $t$ is a palindrome when $t = t_{rev}$ ($t_{rev}$ is a reversed string for $t$).

## More precise statement

It's clear that in the worst case we can have $O(n^2)$ palindrome strings, and at the first glance it seems that there is no linear algorithm for this problem.

But the information about the palindromes can be kept in a more compact way: for each position $i = 0\dots n-1$ we'll find the values $d_1[i]$ and $d_2[i]$, denoting the number of palindromes accordingly with odd and even lengths with centers in the position $i$.

For instance, string $s = abababc$ has three palindromes with odd length with centers in the position $s[3] = b$, i. e. $d_1[3] = 3$:

$a\ \overbrace{b\ a\ \underbrace{b}_{s_3}\ a\ b}^{d_1[3]=3} c$

And string $s = cbaabd$ has two palindromes with even length with centers in the position $s[3] = a$, i. e. $d_2[3] = 2$:

$c\ \overbrace{b\ a\ \underbrace{a}_{s_3}\ b}^{d_2[3]=2} d$

So the idea is that if we have a sub-palindrome with length $l$ with center in some position $i$, we also have sub-palindromes with lengths $l-2$, $l-4$ etc. with centers in $i$. So these two arrays $d_1[i]$ and $d_2[i]$ are enough to keep the information about all the sub-palindromes in the string.

It's a surprising fact that there is an algorithm, which is simple enough, that calculates these "palindromity arrays" $d_1[]$ and $d_2[]$ in linear time. The algorithm is described in this article.

## Solution

In general, this problem has many solutions: with String Hashing it can be solved in $O(n\cdot \log n)$, and with Suffix Trees and fast LCA this problem can be solved in $O(n)$.

But the method described here is sufficiently simpler and has less hidden constant in time and memory complexity. This algorithm was discovered by Glenn K. Manacher in 1975.

## Trivial algorithm

To avoid ambiguities in the further description we denote what "trivial algorithm" is.

It's the algorithm that does the following. For each center position $i$ it tries to increase the answer by one until it's possible, comparing a pair of corresponding characters each time.

Such an algorithm is slow, it can calculate the answer only in $O(n^2)$.

The implementation of the trivial algorithm is:

vector<int> d1(n),  d2(n);
for (int i = 0; i < n; i++) {
d1[i] = 1;
while (0 <= i - d1[i] && i + d1[i] < n && s[i - d1[i]] == s[i + d1[i]]) {
d1[i]++;
}

d2[i] = 0;
while (0 <= i - d2[i] - 1 && i + d2[i] < n && s[i - d2[i] - 1] == s[i + d2[i]]) {
d2[i]++;
}
}

## Manacher's algorithm

We describe the algorithm to find all the sub-palindromes with odd length, i. e. to calculate $d_1[]$. The solution for all the sub-palindromes with even length (i.e. calculating the array $d_2[]$) will be a minor modification for this one.

For fast calculation we'll maintain the borders $(l, r)$ of the rightmost found sub-palindrome (i. e. the palindrome with maximal $r$). Initially we set $l = 0, r = -1$.

So, we want to calculate $d_1[i]$ for the next $i$, and all the previous values in $d_1[]$ have been already calculated. We do the following:

• If $i$ is outside the current sub-palindrome, i. e. $i > r$, we'll just launch the trivial algorithm.

So we'll increase $d_1[i]$ consecutively and check each time if the current rightmost substring $[i - d_1[i]\dots i + d_1[i]]$ is a palindrome. When we find the first mismatch or meet the boundaries of $s$, we'll stop. In this case we've finally calculated $d_1[i]$. After this, we must not forget to update $(l, r)$. $r$ should be updated in such a way that it represents the last index of the current rightmost sub-palindrome.

• Now consider the case when $i \le r$. We'll try to extract some information from the already calculated values in $d_1[]$. So, let's find the "mirror" position of $i$ in the sub-palindrome $(l, r)$, i.e. we'll get the position $j = l + (r - i)$, and we check the value of $d_1[j]$. Because $j$ is the position symmetrical to $i$, we'll almost always can assign $d_1[i] = d_1[j]$. Illustration of this (palindrome around $j$ is actually "copied" into the palindrome around $i$):

$$\ldots\ \overbrace{ s_l\ \ldots\ \underbrace{ s_{j-d_1[j]+1}\ \ldots\ s_j\ \ldots\ s_{j+d_1[j]-1}\ }_\text{palindrome}\ \ldots\ \underbrace{ s_{i-d_1[j]+1}\ \ldots\ s_i\ \ldots\ s_{i+d_1[j]-1}\ }_\text{palindrome}\ \ldots\ s_r\ }^\text{palindrome}\ \ldots$$

But there is a tricky case to be handled correctly: when the "inner" palindrome reaches the borders of the "outer" one, i. e. $j - d_1[j] + 1 \le l$ (or, which is the same, $i + d_1[j] - 1 \ge r$). Because the symmetry outside the "outer" palindrome is not guaranteed, just assigning $d_1[i] = d_1[j]$ will be incorrect: we do not have enough data to state that the palindrome in the position $i$ has the same length.

Actually, we should restrict the length of our palindrome for now, i. e. assign $d_1[i] = r - i + 1$, to handle such situations correctly. After this we'll run the trivial algorithm which will try to increase $d_1[i]$ while it's possible.

Illustration of this case (the palindrome with center $j$ is restricted to fit the "outer" palindrome):

$$\ldots\ \overbrace{ \underbrace{ s_l\ \ldots\ s_j\ \ldots\ s_{j+(j-l)}\ }_\text{palindrome}\ \ldots\ \underbrace{ s_{i-(r-i)}\ \ldots\ s_i\ \ldots\ s_r }_\text{palindrome}\ }^\text{palindrome}\ \underbrace{ \ldots \ldots \ldots \ldots \ldots }_\text{try moving here}$$

It is shown in the illustration that though the palindrome with center $j$ could be larger and go outside the "outer" palindrome, but with $i$ as the center we can use only the part that entirely fits into the "outer" palindrome. But the answer for the position $i$ ($d_1[i]$) can be much bigger than this part, so next we'll run our trivial algorithm that will try to grow it outside our "outer" palindrome, i. e. to the region "try moving here".

Again, we should not forget to update the values $(l, r)$ after calculating each $d_1[i]$.

Also we'll repeat that the algorithm was described to calculate the array for odd palindromes $d_1[]$, the algorithm is similar for the array of even palindromes $d_2[]$. The required modifications can be seen in the code below.

## Complexity of Manacher's algorithm

At the first glance it's not obvious that this algorithm has linear time complexity, because we often run the naive algorithm while searching the answer for a particular position.

However, a more careful analysis shows that the algorithm is linear. In fact, Z-function building algorithm, which looks similar to this algorithm, also works in linear time.

We can notice that every iteration of trivial algorithm increases $r$ by one. Also $r$ cannot be decreased during the algorithm. So, trivial algorithm will make $O(n)$ iterations in total.

Also, other parts of Manacher's algorithm work obviously in linear time. Thus, we get $O(n)$ time complexity.

## Implementation of Manacher's algorithm

For calculating $d_1[]$, we get the following code. Things to note:

• $i$ is the index of the center letter of the current palindrome.
• $d_1[]$ stores the odd palindromes. So, if $i$ exceeds $r$, $k$ is initialized to 1, as a single letter is a palindrome in itself. For $d_2[]$, $k$ will be initialized to 0.
• If $i$ does not exceed $r$, $k$ is either initialized to the $d_1[j]$, where $j$ is the mirror position of $i$ in $(l,r)$, or $k$ is restricted to the size of the "outer" palindrome.
• The while loop denotes the trivial algorithm. We launch it irrespective of the value of $k$.
• If the size of palindrome centered at $i$ is $x$, then $d_1[i]$ stores $(x+1)/2$.
vector<int> d1(n);
for (int i = 0, l = 0, r = -1; i < n; i++) {
int k = (i > r) ? 1 : min(d1[l + r - i], r - i + 1);
while (0 <= i - k && i + k < n && s[i - k] == s[i + k]) {
k++;
}
d1[i] = k--;
if (i + k > r) {
l = i - k;
r = i + k;
}
}

For calculating $d_2[]$, the code looks similar, but with minor changes in arithmetical expressions. Things to note:

• Since an even palidrome will have two centers, $i$ is the latter of the two center indices.
• if $i$ exceeds $r$, $k$ is initialized to 0, as a single letter is not an even palindrome.
• If the size of palindrome centered at $i$ is $x$, then $d_2[i]$ stores $x/2$
vector<int> d2(n);
for (int i = 0, l = 0, r = -1; i < n; i++) {
int k = (i > r) ? 0 : min(d2[l + r - i + 1], r - i + 1);
while (0 <= i - k - 1 && i + k < n && s[i - k - 1] == s[i + k]) {
k++;
}
d2[i] = k--;
if (i + k > r) {
l = i - k - 1;
r = i + k ;
}
}


text:
Generalized Eigenvalue Problem Derivatives
Post
Cancel

# Generalized Eigenvalue Problem Derivatives

The eigenvector problem is ubiquitous in many areas of mathematics, physics and computer science. I recently found myself needing the solution to the generalized eigenvalue problem and discovered an implementation doesn’t exist in jax. While wrapping low-level cuda code is mechanical enough, this doesn’t help with one of the core features of jax: auto-differentiation.

A naive solution would be to re-write the algorithm that generates the decomposition. While this would get us a differentiable solution, it would suffer some major flaws. Most obviously:

• we’d have more code to maintain;
• it likely would not be as efficient as the CUDA implementation; and
• the automatically calculated derivatives would still be sub-optimal.

Instead, we’re going to calculate derivatives using implicit differentiation. This is hardly a new idea, and there are numerous sources detailing results for the standard eigenvalue problem. That said, I was unable to find anything on the generalized problem. In this post, we’ll start by formally defining the problem, before adapting the approach of [Boeddeker et al.][boeddeker] to the generalized problem. In doing so, we identify what we believe to be a flaw in their solution, though identity a work-around for the self-adjoint case. A good understanding of the basics of linear algebra and calculus is assumed.

## Problem Description

We consider the problem of computing partial derivatives of $\W$ and $\V$ given the solution to the generalized eigenvalue and all other relevant partial derivatives. We start with the definition of the solutions to the problem, i.e.

$$\A\W = \B \W \V, \label{eqn:definition}$$ where the columns of $\W$ are the eigenvectors, $\V$ is a diagonal matrix with eigenvalues on the diagonal and square matrices $\A$ and $\B$ define the problem. We limit discussion to the case where the eigenvalues are distinct.

We note there is a degree of freedom in the magnitude of the returned eigenvectors. We ignore this for the moment, but will come back to it.

## A Partial Solution to the General Problem

We begin by taking the partial deriative of Equation \ref{eqn:definition}:

$$\dA \W + \A \dW = \dB \W \V + \B \dW \V + \B \W \dV$$

Pre-multiplying by $\Winv \Binv$ yields $$\Winv \Binv \dA \W + \Winv \Binv \A \dW = \Winv \Binv \left[\B\dW + \dB \W \right] \V + \Winv \Binv \B \W \dV$$

The second term on each side can be simplified: $$\Winv \Binv \dA \W + \V \Winv \dW = \Winv \Binv \left[\B \dW + \dB \W\right] \V + \dV$$

Expanding and simplifying the first term on the right and shifting the second term on the left to the right yields $$\Winv \Binv \dA \W = \Winv \dW \V - \V \Winv \dW + \Winv \Binv \dB \W \V + \dV$$

We can simplify the first two terms on the right by applying the identity $AD - DA = E \circ A$ for a diagonal matrix $D$, where $e_{ij} = d_{ii} - d_{jj}$ and $\circ$ denotes the Hadamard (element-wise) product.

$$\Winv \Binv \dA \W = E \circ \Winv \dW + \Winv \Binv \dB \W \V + \dV \label{eqn:critical}$$

By noting that $I \circ E = \mathbf{0}$, the derivatives of the eigenvalues $\V$ can be decoupled from those of $\W$ by taking the Hadamard product of the identity matrix with Equation \ref{eqn:critical}: $$\dV = I \circ \left[\Winv \Binv \left(\dA \W - \dB \W \V\right)\right].$$

Extracting derivatives of $\W$ is trickier. To start with, we define matrix $F$ of the same size as $E$ with elements given by

$$f_{ij} = \begin{cases} 0 & i = j \\ \frac{1}{e_ij} & i \neq j \end{cases}.$$

By noting the $F \circ E = \mathbf{1} - I$, taking the Hadamard product of equation \ref{eqn:critical} with $F$ yields $$(\mathbf{1} - I)\circ \Winv \dW = F \circ \left[\Winv \Binv \left(\dA \W - \dB \W \V\right)\right]. \label{eqn:troublesome}$$

While everything on the right is straight-forward to compute, the Hadamard product on the left presents difficulties. Giles proposes choosing scales of the eigenvectors such that the diagonal of $\Winv \dW$ are zero, though does not discuss how to compute such a scale. Note that contrary to the claims of Boeddeker et al., using normalized eigenvectors (or more generally, eigenvectors of fixed magnitude) does not guarantee this property.

We make no effort to resolve this for the general case. We argue that use-cases of eigenvectors should be invariant to this scaling, and derivatives would likely be computed more easily with respect to these quantities instead. For the interested reader, a good discussion of the problem is given here.

One special case where computing derivatives is with respect to eigenvectors is relatively straight forward is in the self-adjoint case, where $\A = \A^H$ and $\B = \B^H$. In this case, eigenvectors form an orthogonal basis. Equation \ref{eqn:troublesome} can then be computed by enforcing magnitudes such that the eigenvectors are orthogonal with respect to $\B$, i.e. $$\W^H \B \W = I \iff{\Winv = \W^H \B}. \label{eqn:orthonormal}$$

This allows for some simplifications to Equation \ref{eqn:troublesome}: $$(\mathbf{1} - I)\circ \W^H \B \dW = F \circ \left[\W^H \left(\dA \W - \dB \W \V\right)\right]. \label{eqn:less-troublesome}$$

We still have the problem of the identity matrix in the Hadamard operator. We can distribute the product operator over addition ($[\mathbf{1} + I] \circ X = X + I \circ X$) but this still isn’t in a form we can compute easily compute. Fortunately, the orthonormal constraint gives us a mechanism by which we can express the diagonal values of $\dW$ in terms of other known values.

To understand how, we begin by differentiating Equation \ref{eqn:orthonormal} $$\dWH \B \W + \W^H \dB \W + \W^H \B \dW = \mathbf{0}.$$

Since the first and third term are Hermitian transposes of each other, their complex components will cancel out on the diagonal. This means we can expression the real part of our troublesome term in terms of less problematic components.

$$I \circ \Re\left[\W^H \B \dW\right] = -\frac{1}{2} I \circ \W^H \dB \W. \label{eqn:identity-hadamard}$$

Alas, this does not help us in evaluating the imaginary part. We do still have some freedom however: it can be shown if $\V, \W$ is a solution to the problem and satisfies the described orthogonality constraint, then $\V \mathbf{R}(\Theta)$ is also a solution, where $\mathbf{R}(\Theta)$ is a diagonal matrix with entries given by rotations in the complex plan, $R_{jj} = e^{i\theta_j}$. Exactly how we define $\theta_j$ in such a way that makes the remaining troublesome term computable is left to the reader (i.e. I’ve spent entirely too long trying and got nowhere).

## Real Symmetric Case

Like we did above, rather than overcome this dilemna we’ll simply claim we were never interested in the solution to the complex problem anyway.

For real $\A$ and $\B$ the imaginary part is zero. Substituting Equation \ref{eqn:identity-hadamard} into Equation \ref{eqn:less-troublesome} and rearranging gives us our final expression for the eigenvector derivatives: $$\dW = \W \left( F \circ \left[\W^H\left(\dA \W - \dB \W \V\right)\right] - \frac{1}{2} I \circ \left[\W^H \dB \W\right]\right).$$

## Jax Implementation

No amount of reviewing maths will ever make me trust a solution as much as numerical validation. We implemented the above in jax (also available in this gist).

eigh_impl.py:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 """Versions based on 4.60 and 4.63 of https://arxiv.org/pdf/1701.00392.pdf .""" import jax import jax.numpy as jnp import numpy as np def _T(x): return jnp.swapaxes(x, -1, -2) def _H(x): return jnp.conj(_T(x)) def symmetrize(x): return (x + _H(x)) / 2 def standardize_angle(w, b): if jnp.isrealobj(w): return w * jnp.sign(w[0, :]) else: # scipy does this: makes imag(b[0] @ w) = 1 assert not jnp.isrealobj(b) bw = b[0] @ w factor = bw / jnp.abs(bw) w = w / factor[None, :] sign = jnp.sign(w.real[0]) w = w * sign return w @jax.custom_jvp # jax.scipy.linalg.eigh doesn't support general problem i.e. b not None def eigh(a, b): """ Compute the solution to the symmetrized generalized eigenvalue problem. a_s @ w = b_s @ w @ np.diag(v) where a_s = (a + a.H) / 2, b_s = (b + b.H) / 2 are the symmetrized versions of the inputs and H is the Hermitian (conjugate transpose) operator. For self-adjoint inputs the solution should be consistent with scipy.linalg.eigh i.e. v, w = eigh(a, b) v_sp, w_sp = scipy.linalg.eigh(a, b) np.testing.assert_allclose(v, v_sp) np.testing.assert_allclose(w, standardize_angle(w_sp)) Note this currently uses jax.linalg.eig(jax.linalg.solve(b, a)), which will be slow because there is no GPU implementation of eig and it's just a generally inefficient way of doing it. Future implementations should wrap cuda primitives. This implementation is provided primarily as a means to test eigh_jvp_rule. Args: a: [n, n] float self-adjoint matrix (i.e. conj(transpose(a)) == a) b: [n, n] float self-adjoint matrix (i.e. conj(transpose(b)) == b) Returns: v: eigenvalues of the generalized problem in ascending order. w: eigenvectors of the generalized problem, normalized such that w.H @ b @ w = I. """ a = symmetrize(a) b = symmetrize(b) b_inv_a = jax.scipy.linalg.cho_solve(jax.scipy.linalg.cho_factor(b), a) v, w = jax.jit(jax.numpy.linalg.eig, backend="cpu")(b_inv_a) v = v.real # with loops.Scope() as s: # for _ in s.cond_range(jnp.isrealobj) if jnp.isrealobj(a) and jnp.isrealobj(b): w = w.real # reorder as ascending in w order = jnp.argsort(v) v = v.take(order, axis=0) w = w.take(order, axis=1) # renormalize so v.H @ b @ H == 1 norm2 = jax.vmap(lambda wi: (wi.conj() @ b @ wi).real, in_axes=1)(w) norm = jnp.sqrt(norm2) w = w / norm w = standardize_angle(w, b) return v, w @eigh.defjvp def eigh_jvp_rule(primals, tangents): """ Derivation based on Boedekker et al. https://arxiv.org/pdf/1701.00392.pdf Note diagonal entries of Winv dW/dt != 0 as they claim. """ a, b = primals da, db = tangents if not all(jnp.isrealobj(x) for x in (a, b, da, db)): raise NotImplementedError("jvp only implemented for real inputs.") da = symmetrize(da) db = symmetrize(db) v, w = eigh(a, b) # compute only the diagonal entries dv = jax.vmap( lambda vi, wi: -wi.conj() @ db @ wi * vi + wi.conj() @ da @ wi, in_axes=(0, 1), )(v, w) dv = dv.real E = v[jnp.newaxis, :] - v[:, jnp.newaxis] # diagonal entries: compute as column then put into diagonals diags = jnp.diag(-0.5 * jax.vmap(lambda wi: wi.conj() @ db @ wi, in_axes=1)(w)) # off-diagonals: there will be NANs on the diagonal, but these aren't used off_diags = jnp.reciprocal(E) * (_H(w) @ (da @ w - db @ w * v[jnp.newaxis, :])) dw = w @ jnp.where(jnp.eye(a.shape[0], dtype=np.bool), diags, off_diags) return (v, w), (dv, dw)

Dirty testing script:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 from eigh_impl import symmetrize, eigh, standardize_angle import jax.numpy as jnp import numpy as np import jax.test_util as jtu import scipy.linalg jnp.set_printoptions(3) rng = np.random.default_rng(0) n = 5 is_complex = False def make_spd(x): n = x.shape[0] return symmetrize(x) + n * jnp.eye(n) def get_random_square(rng, size, is_complex=True): real = rng.uniform(size=size).astype(np.float32) if is_complex: return real + rng.uniform(size=size).astype(np.float32) * 1j return real a = make_spd(get_random_square(rng, (n, n), is_complex=is_complex)) b = make_spd(get_random_square(rng, (n, n), is_complex=is_complex)) vals, vecs = eigh(a, b) # ensure solution satisfies the problem np.testing.assert_allclose(a @ vecs, b @ vecs @ jnp.diag(vals), atol=1e-5) # ensure vectors are orthogonal w.r.t b np.testing.assert_allclose(vecs.T.conj() @ b @ vecs, jnp.eye(n), atol=1e-5, rtol=1e-5) # ensure eigenvalues are ascending np.testing.assert_array_less(vals[:-1], vals[1:]) jtu.check_grads(eigh, (a, b), 2, modes=["fwd"]) # ensure values consistent with scipy vals_sp, vecs_sp = scipy.linalg.eigh(a, b) print("scipy") print(vecs_sp) print("this work") print(vecs) np.testing.assert_allclose(vals, vals_sp, rtol=1e-4, atol=1e-5) np.testing.assert_allclose(vecs, standardize_angle(vecs_sp, b), rtol=1e-4, atol=1e-5) print("success")


text:
# Étude de la quadrangulation infinie uniforme

Abstract : Quadrangulations are proper embeddings of finite connected graphs in the two-dimensional sphere for which all faces have degree 4. Our main object of study is the so-called uniform infinite quadrangulation. This infinite random map has been defined in two different ways. The first definition consists in taking the local limit of large random (finite) quadrangulations whose laws are uniform over the set of all quadrangulations with a given size. The second definition takes advantage of a bijection between quadrangulations and well-labelled trees. The starting point is a suitably defined uniform infinite random tree whose law is mapped to the set of all infinite quadrangulations using an extended version of the bijection. In the second chapter of this manuscript, we prove that the two previous definitions lead to the same object. This fact is not trivial as, in the infinite setting, the bijection between quadrangulations and trees is not continuous for the topology of local convergence. The proof depends on studying some combinatorial properties of the bijection and estimating the distribution of vertices with small labels in high generations of the uniform infinite tree. In chapter 3, we use the equivalence of the two definitions to compute scaling limits for the uniform infinite quadrangulation. Indeed, quantities such as the volume of balls around a distinguished point of the quadrangulation can be evaluated by studying the uniform infinite tree. The main technical ingredient of this chapter is the convergence of the rescaled contour functions of the uniform infinite tree towards a stochastic process linked to the Brownian snake.
Mots-clés :
Document type :
Theses
Domain :
Complete list of metadatas

https://tel.archives-ouvertes.fr/tel-00467174
Contributor : Laurent Ménard <>
Submitted on : Friday, March 26, 2010 - 10:49:57 AM
Last modification on : Wednesday, March 27, 2019 - 1:28:10 AM
Long-term archiving on : Friday, October 19, 2012 - 10:40:38 AM

### Identifiers

• HAL Id : tel-00467174, version 1

### Citation

Laurent Ménard. Étude de la quadrangulation infinie uniforme. Mathématiques [math]. Université Pierre et Marie Curie - Paris VI, 2009. Français. ⟨tel-00467174⟩

Record views

Files downloads


text:
1. ## Curves

I'm not sure what the curve of x=-(sqrt(1-(y-1)^2)) would look like....anyone know??

2. Originally Posted by s7b
I'm not sure what the curve of x=-(sqrt(1-(y-1)^2)) would look like....anyone know??
Square both sides and add (y-1)^2 to both sides and you should see a familiar looking form.

3. oh ok, so i think it would be a circle of radius one centered about (0,1) ??

4. Originally Posted by s7b
oh ok, so i think it would be a circle of radius one centered about (0,1) ??
Sounds good. But you weren't given the ending formula, you were given the negative square root part. How does that affect the graph?

5. Would that mean that it is reflected in the x-axis giving a circle of radius one centered at (0,-1)

6. When you take conics and apply a square root, the result is a split in the graph. In your problem x is always going to be negative so I believe this implies that it is the circle you first mentioned, cut in half vertically, and the left side of the half circle.