999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

本質(zhì)擬對稱非負(fù)張量H-譜半徑的算法

2025-03-21 00:00:00林志興呂洪斌

摘要: 首先, 定義一類本質(zhì)擬對稱非負(fù)張量, 其是涵蓋本質(zhì)正張量、 弱正張量和廣義弱正張量的更廣泛的張量類. 其次, 應(yīng)用張量的H-特征值在對角相似變換下的不變性質(zhì), 給出本質(zhì)擬對稱非負(fù)張量H-譜半徑的算法, 并用數(shù)值例子說明算法的有效性.

關(guān)鍵詞: 本質(zhì)擬對稱非負(fù)張量; H-譜半徑; 算法

中圖分類號: O183.2文獻(xiàn)標(biāo)志碼: A文章編號: 1671-5489(2025)02-0411-06

Algorithm for H-Spectral Radius of EssentiallyQuasi-symmetric Non-negative Tensors

LIN Zhixing, L Hongbin

(Fujian Key Laboratory of Financial Information Processing (Putian University), Putian 351100, Fujian Province, China)

Abstract: Firstly, we defined a class of essentially quasi-symmetric non-negative tensors, which encompassed a broader class of tensors covering essentially positive

tensors, weakly positive tensors and generalized weakly positive tensors. Secondly, we gave an algorithm for the H-spectral radius of an essentially quasi-symmetric non-negati

ve tensorby applying the property that the H-eigenvalues of the tensor were invariant under diagonal similarity transformations, and usednumerical examples

toillustrate the effectiveness of the algorithm.

Keywords: essentially quasi-symmetric non-negative tensor; H-spectral radius; algorithm

收稿日期: 2024-07-12.

第一作者簡介: 林志興(1973—), 男, 漢族, 碩士, 教授, 從事矩陣?yán)碚摷捌鋺?yīng)用的研究, E-mail: linzhixing@163.com.通信作者簡介: 呂洪斌

(1964—), 男, 漢族, 博士, 教授, 從事數(shù)值代數(shù)、 矩陣?yán)碚摷捌鋺?yīng)用的研究, E-mail: hbinlv@126.com.

基金項(xiàng)目: 國家自然科學(xué)基金(批準(zhǔn)號: 61772292)、 福建省自然科學(xué)基金(批準(zhǔn)號: 2023J01997)和莆田市科技項(xiàng)目(批準(zhǔn)號: 2023SZ3001PTXY17).

張量是矩陣的一種高維推廣, 在信號和圖像處理、 連續(xù)介質(zhì)物理、 數(shù)據(jù)挖掘和處理、 非線性優(yōu)化、 物理學(xué)的彈性分析和高階統(tǒng)計(jì)學(xué)等領(lǐng)域應(yīng)用廣泛1-3].一個m階n維實(shí)張量是nm個元素的多維數(shù)組

A=(ai1i2…im), ai1i2…im∈瘙綆, ij=1,2,…,n, j=1,2,…,m.

用瘙綆[m,n表示實(shí)數(shù)域瘙綆上的m階n維張量的集合. 若ai1i2…im≥0(ij=1,2,…,n, j=1,2,…,m), 則稱A為非負(fù)

張量, 記為A∈瘙綆[m,n+. 此外, 本文中瘙綇n×n和瘙綆n×n+分別表示復(fù)數(shù)域和實(shí)數(shù)域上的n×n階矩陣及其非負(fù)矩陣集合,

瘙綆n+和瘙綆n++分別表示n維歐氏空間上的非負(fù)非零向量和正向量集合.

目前, 關(guān)于張量特征值的研究已有很多結(jié)果: Qi[4和Lim[5分別定義了張量的特征值; 文獻(xiàn)[6]給出了不可約非負(fù)張量H-譜半徑的NQZ算法;文獻(xiàn)[7]證明了算法對本質(zhì)正張量的收斂性; 文獻(xiàn)[8]提出了非負(fù)張量H-譜半徑的LZI算法, 并證明了LZI算法對一類不可約非負(fù)張量的收斂性; 文獻(xiàn)[9]證明了弱正張量H-譜半徑NQZ算法的線性收斂性; 文獻(xiàn)[10]給出了廣義弱正張量的定義和H-譜半徑算法的線性收斂性; 文獻(xiàn)[11]給出了弱本質(zhì)不可約非負(fù)張量H-譜半徑的擬冪型算法. 本文定義一類本質(zhì)擬對稱非負(fù)張量, 并給出不可約本質(zhì)擬對稱非負(fù)張量H-譜半徑的對角相似算法.

1 相關(guān)定義與結(jié)果

定義1[12設(shè)A=(ai1i2…im)∈瘙綆[m,n, 若有非空真子集J〈n〉={1,2,…,n}, 使得aii2…im=

0, i∈J, i2,…,im∈〈n〉\J, 則稱張量A是可約的, 否則稱其為不可約的.

定義2[7設(shè)A=(ai1i2…im)∈瘙綆[m,n+, 則稱M(A)=(aij)

∈瘙綆n×n+是張量A的優(yōu)化矩陣, 其中aij=aij…j, i,j∈〈n〉.

由定義1和定義2知, 若張量A=(ai1i2…im)∈瘙綆[m,n+的優(yōu)化矩陣M(A)是不可約的, 則A是不可約張量.

定義3[7,9-10設(shè)A=(ai1i2…im)∈瘙綆[m,n+, M(A)=(aij

)∈瘙綆n×n+是張量A的優(yōu)化矩陣. 若aijgt;0(i,j∈〈n〉), 則稱A是本質(zhì)正張量; 若aijgt;0

(i,j∈〈n〉, i≠j), 則稱A是弱正張量; 若存在i0∈〈n〉, 使得ai0jgt;0, aji0gt;0(i0,j∈〈n〉, j≠i0), 則稱A是廣義弱正張量.

本文推廣本質(zhì)正張量、 弱正張量和廣義弱正張量類, 定義一類本質(zhì)擬對稱非負(fù)張量, 并給出其H-特征值的算法.

定義4 設(shè)A=(ai1i2…im)∈瘙綆[m,n+, M(A)=(aij)∈瘙綆

n+是張量A的優(yōu)化矩陣. 若aijgt;0.ajigt;0(i≠j), 則稱A是本質(zhì)擬對稱非負(fù)張量;

圖1 非負(fù)張量類的包含關(guān)系Fig.1 Inclusion relations of non-negative tensor classes

若M(A)是不可約的, 則稱A是不可約本質(zhì)擬對稱非負(fù)張量.

由定義4知, 本質(zhì)正張量、 弱正張量和廣義弱正張量均是特殊的本質(zhì)擬對稱非負(fù)張量類, 其與本原張量和不可約非負(fù)張量之間的包含關(guān)系如圖1所示.

定義5[4-5設(shè)A=(ai1i2…im)

∈瘙綆[m,n, 若存在λ∈瘙綇和非零向量x=(x1,x2,…,xn)T, 使得

Axm-1=λx[m-1

則λ稱為張量A的特征值, x為張量A的相應(yīng)于特征值λ的特征向

量, 其中Axm-1和x[m-1為n維向量, 它們的第i個元素分別為

(Axm-1)i=∑ni2,…,im=1a

ii2…imxi2…xim,(x[m-1)i=xm-1i.

若λ和x的元素都是實(shí)的, 則λ稱為張量A的H-特征值, x稱為相應(yīng)于

特征值λ的H-特征向量. 本文用σH(A)表示張量A的H-特征

值集合, ρ(A)=maxλ∈σH(A)λ稱為張量A的H-譜半徑.

定義6[13設(shè)A=(ai1i2…im)∈瘙綆[m,n, B=(bi1i2…im)∈瘙綆

[m,n, 如果存在對角可逆矩陣D=diag(d1,d2,…,dn)∈瘙綆n×n, 使得B=D-(m-1)·A·

D, 則稱張量A,B是對角相似的, 其中bi1i

2…im=ai1i2…imd-(m-1)i1di2…dim, ij=1,2,…,n, j=1,2,…,m.

定理1[13設(shè)A=(ai1i2…im)∈瘙綆[m,n, B=(bi1i2…i

m)∈瘙綆[m,n, 如果張量 A,B是對角相似的, 則σH(A)=σH(B).

類似于非負(fù)矩陣, 非負(fù)張量的H-譜半徑有如下基本性質(zhì).

定理2[12設(shè)A=(ai1i2…im)∈瘙綆[m,n+, 則:

1) ρ(A)是張量A的特征值, 且ρ(A)有相應(yīng)的非負(fù)特征向量

x∈瘙綆n+, 如果不計(jì)倍數(shù), 則x∈瘙綆n+是唯一的;

2) 若張量A是不可約的, 則ρ(A)有相應(yīng)的正特征向量x∈瘙綆n++.

Chang等12將關(guān)于非負(fù)矩陣譜半徑上下界的經(jīng)典Perron-Frobenius定理推廣到非負(fù)張量上, 得到如下結(jié)果:

定理3[12設(shè)張量A=(ai1i2…im)∈瘙綆[m,n+, 則

mini∈〈n〉 ri(A)≤ρ(A)≤maxi∈〈n〉 ri(A),

其中ri(A)=∑ni2,…,im=1aii2…im, i∈〈n〉.

設(shè)矩陣A=(aij)∈瘙綇n×n, 記A的有向圖14為Γ(A).

如果i1,i2,…,ir+1互不相同, 且aisis+1≠0(s=1,2,…,r), 則稱ei1→ir+1是Γ(A)上的有向路徑, 并稱r是有向路徑ei1→ir+1的長度.

定義7[14設(shè)矩陣A=(aij)∈瘙綇n×n, Γ(A)是A的有向圖. 如

果對任意兩個不同的結(jié)點(diǎn)i,j∈〈n〉都有i到j(luò)的有向路徑, 則稱有向圖Γ(A)是強(qiáng)連通的.

定理4[14設(shè)矩陣A=(aij)∈瘙綇n×n, 則A是不可約的當(dāng)且僅當(dāng)A的有向圖Γ(A)是強(qiáng)連通的.

2 算法構(gòu)造

設(shè)張量A=(ai1i2…im)∈瘙綆[m,n+, 取-mini∈〈n〉 ai…ilt;θlt;maxi∈〈n〉 ri(A), 記A=∶A(0)=(a(0)i1i2…im)ni1,i2,…,im=1. 對k=0,1,2,…, 記

r(k)i=∑i2,…,im∈〈n〉a(k)ii2…im, r(k)=maxi∈〈n〉

r(k)i, r(k)=mini∈〈n〉 r(k)i,Dk=1(r(k)+

θ)1/(m-1)[diag(r(k)1,r(k)2,…,r(k)n)+θI]1/(m-1)

A(k+1)=(Dk)-(m-1)·A(k)·Dk,

其中I為單位矩陣. 則有非負(fù)張量H-譜半徑的如下算法.

算法1

設(shè)張量A=(ai1i2…im)∈瘙綆[m,n+, εgt;0, k=0.

步驟1) 計(jì)算r(k)i=∑i2,…,im∈〈n〉a(k)ii2…im, r(k)=maxi∈〈n〉 r(k)

i, r(k)=mini∈〈n〉 r(k)i.

步驟2) 如果r(k)-r(k)lt;ε, 轉(zhuǎn)步驟3); 否則, 記

Dk=1(r(k)+θ)1/(m-1)[diag(r(k)1,r(k)2,…,r(k)n)+θI]1/(m-1),A

(k+1)=(Dk)-(m-1)·A(k)·Dk.

令k∶=k+1, 轉(zhuǎn)步驟1).

步驟3) 輸出ρ(A)=12(r(k)+r(k)).

3 算法的收斂性

下面證明算法1對不可約本質(zhì)擬對稱非負(fù)張量H-譜半徑的計(jì)算是收斂的.

引理1 設(shè)張量A=(ai1i2…im)∈瘙綆[m,n+, 則r(k)單調(diào)遞增有上界, r(k)單調(diào)遞減有下界.

證明: 應(yīng)用定理3和定理1, 對任意的k∈瘙綄+∪{0}, 有

r(0)≤r(k+1)=maxi∈〈n〉 r(k+1)i=maxi∈〈

n〉∑i2,…,im∈〈n〉(a(k)ii2…im+θ)(r(k)i2+θ)1/(m-1)·…·(r

(k)im+θ)1/(m-1)r(k)i+θ-θ≤(r(k)+θ)∑i2,…,i

m∈〈n〉(a(k)ii2…im+θ)r(k)i+θ-θ=r(k)

因此, r(k)單調(diào)遞減有下界. 類似可證r(k)單調(diào)遞增有上界. 證畢.

引理2 設(shè)A=(ai1i2…im)∈瘙綆[m,n+是不可約本質(zhì)擬對稱張量, 則對任意的k∈瘙綄+和任意的a(0)

ij…jgt;0(i≠j), 有

a(k)ij=a(k)ij…j≥a∶=a2/r(0)gt;0,

其中a=min{aij…jgt;0: i,j∈〈n〉, i≠j}.

證明: 顯然, 對任意的k∈瘙綄+, 若aij…jgt;0, 則a(k)ij…jgt;0. 由引理1有

r(0)≥r(k)≥a(k)ij…j=a(k-1)ij…j·r(k-1)j+θr(k-1)i+θ=…=a(0)ij…j·∏k-1t=0

(r(t)j+θ)∏k-1t=0(r(t)i+θ)≥a·∏k-1

t=0(r(t)j+θ)∏k-1t=0(r(t)i+θ),(1)

從而有∏k-1t=0(r(t)j+θ)∏k-1t=0(r(t)i+θ)≤r(0)a.

又由∏k-1t=0(r(t)j+θ)∏k-1t=0(r(t)i+θ)·∏k-1t=0(r(t)i+θ)∏k-1t=0(r(t)j+θ)=1

及aji…igt;0, 有

∏k-1t=0(r(t)j+θ)∏k-1t=0(r(t)i+θ)≥a2/r(0)

再由式(1)有a(k)ij=a(k)ij…j≥a∶=a2/r(0)gt;0. 證畢.

定理5 設(shè)A=(ai1i2…im)∈瘙綆[m,n+是不可約本質(zhì)擬對稱張量, 則對k∈瘙綄+∪{0}, 由算法1有

r((k+1)n)-r((k+1)n)≤α(r(kn)-r(kn)),

其中0lt;α=1-a0r(0)+θnlt;1, a0∶=mini∈〈n〉{a,aii…i+θ}gt;

0, 進(jìn)而有 limk→∞ r(k)=limk→∞ r(k) =ρ(A).

證明: 設(shè)r(0)lt;r(0), 否則, 由定理3知r(0)=r(0)=ρ(A). 由

張量序列{A(k)}∞k=0的構(gòu)造過程知, 當(dāng)a(0)i1i2…im≠0時, a(k)i1i2…im≠0, i1,i2

,…,im∈〈n〉; a(k)i…i=a(0)i…i=ai…i, i∈〈n〉.

設(shè)A=(ai1i2…im)∈瘙綆[m,n+是不可約本質(zhì)擬對稱張量, 定義A(k)+θE=∶(k)=(a(k)i1i2…im), k=0,1,2,…,

其中E∈瘙綆[m,n+是單位張量. 則由引理2知, 對任意的a(0)ij…jgt;0, i,j∈〈n

〉, k∈瘙綄+, 有a(k)ij…j≥a0gt;0.

取i0∈E0={i∈〈n〉: ri(A)=mini∈〈n〉 ri(A)}, 則有

r(1)i0=∑ni2,…,im=1a(1)i0i2…im=r(0)-∑ni2,…,i

m=1a(0)i0i2…im(r(0)+θ)-∏ml=2(r(0)il+θ)

1/(m-1)r(0)i0+θ≤r(0)-a(0)i0…i0(r

(0)+θ)-(r(0)i0+θ)r(0)i0+θ≤r(0)-a0r(0)

+θ·(r(0)-r(0)).(2)

類似地, 利用式(2)并依次推導(dǎo)可得

r(t)i0≤r(0)-a0r(0)+θt(r(0)-r(0)),t=1,2,…,n.(3)

對任意的i∈〈n〉\E0, 由A是不可約本質(zhì)擬對稱張量知, Γ(M(A))是強(qiáng)

連通的, 故存在a1i0…i0gt;0, a21…1gt;0, …, ar-1

r-2…r-2gt;0, air-1…r-1gt;0, 其中i0,1,…,r-1,i互不相同, 且r≤n-1. 由式(3)有

r(2)1≤r(0)-a(1)1i0…i0(r(0)+θ)-(r(1)i0+θ)r(1)1+θ≤r(0)-a0r

(0)+θ·(r(0)-r(1)i0)lt;r(0)-a0r

(0)+θ2(r(0)-r(0)).(4)

類似地, 利用式(4)并依次推導(dǎo)可得

r(t)t-1≤r(0)-a(t-1)t-1t-2…t-2(r(0)+θ)-(r(t-

1)t-2+θ)r(t-1)t-1+θ≤r(0)-a0r

(0)+θ·(r(0)-r(t-1)t-2)lt;r(0)-a0r(0)+θt(r(0)-r(0)), t=3,4,…,n.

從而有r(n)i≤r(0)-a0r(0)+θn(r(0)-r(0)).

于是有

r(n)-r(n)≤maxi∈〈n〉 r(n)i-r(0)

1-a0r(0)+θn(r(0)-r(0)).

類似上述證明, 對任意的k∈瘙綄+, 有

r(kn)-r(kn)≤1-a0r(0)+θn(r((k-1)n)-r

((k-1)n))≤…≤1-a0r(0)+θnk(r(0)-r(0)).

由于0lt;α∶=1-a0r(0)+θnlt;1, 應(yīng)用引理1并令k→∞, 則有

limk→∞ r(kn)-limk→∞ r(kn)=

limk→∞(r(kn)-r(kn))=limk→∞(αk(r(0)-r(0)))=0,

即ρ(A)=limk→∞ r(kn)=limk→∞ r(kn). 證畢.

注1 算法1中引入了參數(shù)θ, 相當(dāng)于一個移位算法, 在設(shè)計(jì)上是有意義的. 同時, 參數(shù)θ在其取值范圍內(nèi)對算法的計(jì)算效率有一定影響, 但還不能從理論上說明θ的最優(yōu)取值.

注2 定理5的條件是非負(fù)張量A是不可約本質(zhì)擬對稱的, 但條件過于嚴(yán)

格. 事實(shí)上, 若記M(A)=(aij)∈瘙綆n×n+, 只要有矩陣(A)=(aij)∈瘙綆n×n+,

使得矩陣(A)擬對稱且不可約, 則由定理5的證明即知算法1收斂, 其中a

ij= aij或0,aij≠0, 0,aij=0.

例如, 對于張量A=(aijk)∈瘙綆[3,4+, 設(shè)M(A)=0

100302201034010, 其余元素任意非負(fù), 取(A)=0

100302001030010, 則

(A)是擬對稱的, 此時算法1收斂.

注3 設(shè)A=(ai1i2…im)∈瘙綆[m,n+是不可約本質(zhì)擬對稱張量, 則由算法1知, x=∏k→∞D(zhuǎn)(k)e即為ρ(A)對應(yīng)的正特征向量, 其中e=(1,1,…,1)T∈瘙綆n++.

4 數(shù)值算例

下面用實(shí)例比較算法1和LZI算法8的計(jì)算效率.

例1 設(shè)張量A=(aijk)∈瘙綆[3,n+, 其中aijj=(i+j)/(2i), j=i+1, i=j+1, i,j=1,2,…,n-1, 其余元素為零, 則有

M(A)=032000340540005600002n-12(n-1)0002n-12n0,

因此, A是不可約本質(zhì)擬對稱非負(fù)張量. 用算法1(取θ=(r(0)-r(0))/2=0.25)和LZI算法分別計(jì)算ρ(A)的迭代步數(shù)和CPU時間, 比較結(jié)果列于表1.

由表1可見, 計(jì)算例1的非負(fù)張量H-譜半徑的算法1 迭代次數(shù)明顯少于LZI算法, 但在CPU用時上有些不足.

綜上所述, 本文定義的本質(zhì)擬對稱非負(fù)張量類, 是本質(zhì)正張量、 弱正張量和廣義弱正張量類的有益推廣, 給出的H-譜半徑的帶有參數(shù)變換的算法, 解決了某些非負(fù)張量H-譜半徑NQZ算法

[6]不收斂、 LZI算法8不能直接計(jì)算的問題, 說明了參數(shù)選擇對算法計(jì)算效率的影響, 并比較了其與LZI算法的效率, 數(shù)值計(jì)算結(jié)果表明本文算法是有意義的.

參考文獻(xiàn)

[1]CICHOCKI A, MANDIC D P, PHAN A H, et al. Tensor Decompositions fo

r Signal Processing Applications: From Two-Way to Multiway Component Analysis[J].IEEE Signal Processing Magazine, 2015, 32(2): 145-163.

[2]SCHULTZ T, SEIDEL H P. Estimating Crossing Fibers: A Tensor Decomposition Ap

proach[J].IEEE Transactions on Visualization and Computer Graphics, 2008, 14(6): 1635-1642.

[3]QI L Q, WANG Y J, WU E X. D-Eigenvalues of Diffusion Kurtos

is Tensors[J].Journal of Computational and Applied Mathematics, 2008, 221(1): 150-157.

[4]QI L Q. Eigenvalues of a Real Supersymmetric Tensor[J].Journal of Symbolic Computation, 2005, 40(6): 1302-1324.

[5]LIM L H. Singular Values and Eigenvalues of Tensors: A Variational Approach[C]

//1st IEEE International Workshop on Computational Advances in Multi-sensor Adaptive Processing. Puerto Vallarta, Mexico:[s.n.], 2005: 129-132.

[6]NG M, QI L Q, ZHOU G L. Finding the Largest Eigenval

ue of a Nonnegative Tensor[J].SIAM Journal of Matrix Analysis and Applications, 2009, 31(3): 1090-1099.

[7]PEARSON K J. Essentially Positive Tensors[J].International Journal of Algebra, 2010, 4(9/10/11/12): 421-427.

[8]LIU Y J, ZHOU G L, IBRAHIM N F. An Always Convergent Algorithm for the Largest

Eigenvalue of an Irreducible Nonnegative Tensor[J].Journal of Computational and Applied Mathematics, 2010, 235(1): 286-292.

[9]ZHANG L P, QI L Q, XU Y. Linear Convergence of the LZI Algorithm for Weakly Positive Tensors [J]. Journal of Computational Mathematics,

2012, 30(1): 24-33.

[10]ZHANG J L, BU C J. An Iterative Method for Finding the Spectral Radius of an I

rreducible Nonnegative Tensor[J].Computational amp; Applied Mathematics, 2021, 40(1): 8-1-8-14.

[11]LIU G M, L H B. An Algorithm for the Spectral Radius of Weakly Essent

ially Irreducible Nonnegative Tensors[J].Calcolo, 2024, 61(1): 8-1-8-19.

[12]CHANG K C, PEARSON K, ZHANG T. Perron-Frobenius Theor

em for Nonnegative Tensors[J].Communications in Mathematical Sciences, 2008, 6(2): 507-520.

[13]SHAO J Y. A General Product of Tensors with Applications[J].Linear Algebra and Its Applications, 2013, 439(8): 2350-2366.

[14]HORD R A, JOHNSON C R. Matrix Analysis[M].2nd ed. Cambridge: Cambridge University Press, 2014: 399-400.

(責(zé)任編輯: 趙立芹)

主站蜘蛛池模板: 欧美日韩中文字幕在线| 亚洲一级毛片免费看| 欧美a在线| 99视频精品全国免费品| 毛片免费网址| 天天躁狠狠躁| 91亚洲国产视频| 亚洲第一成年网| 伊人大杳蕉中文无码| www.精品国产| 原味小视频在线www国产| 在线a网站| AV天堂资源福利在线观看| 国产午夜精品一区二区三| 亚洲婷婷丁香| 国产国拍精品视频免费看| 在线观看亚洲国产| 欧美一级色视频| 国产男人的天堂| 亚洲精品国产乱码不卡| 亚洲激情区| 欧美人与动牲交a欧美精品| 午夜视频免费一区二区在线看| 亚洲欧美成人在线视频| 欧美午夜网| 亚洲欧美精品一中文字幕| 99九九成人免费视频精品| 日韩成人午夜| 亚洲无码视频喷水| 亚洲系列无码专区偷窥无码| 青青青视频91在线 | 免费观看精品视频999| 国产精品成人啪精品视频| 91视频国产高清| 97视频在线观看免费视频| 国产精品嫩草影院视频| 欧美精品不卡| 免费视频在线2021入口| 久草网视频在线| 久久黄色小视频| 女人18毛片一级毛片在线| 重口调教一区二区视频| 东京热av无码电影一区二区| 久久这里只有精品免费| 国产成人啪视频一区二区三区 | 欧美日韩v| 伊人色婷婷| 国产成本人片免费a∨短片| 日韩在线播放欧美字幕| 亚洲女同欧美在线| 五月婷婷综合网| 精品国产成人高清在线| 久久综合结合久久狠狠狠97色| 在线国产91| 九九九国产| 国产91av在线| 高潮爽到爆的喷水女主播视频| 91探花在线观看国产最新| 亚洲国产一区在线观看| 91极品美女高潮叫床在线观看| 午夜国产大片免费观看| 女同国产精品一区二区| 香蕉久久国产超碰青草| 全部免费特黄特色大片视频| 欧美黑人欧美精品刺激| 成年片色大黄全免费网站久久| 91午夜福利在线观看精品| 国产人人乐人人爱| 在线五月婷婷| 99精品国产高清一区二区| 国产青青草视频| 国产视频只有无码精品| 日本道综合一本久久久88| 2018日日摸夜夜添狠狠躁| 国产成人在线小视频| 国产96在线 | 色婷婷综合激情视频免费看| 成人午夜久久| 欧美啪啪视频免码| 欧美性精品不卡在线观看| 久久久久国产精品免费免费不卡| 免费在线色|