机器学习hierarchical clustering_材料学+AI：非监督学习预测新型固态锂离子导体材料...

2019年11月20日，丰田北美研究院的凌晨博士与马里兰大学莫一非教授合作，使用非监督学习的方法，预测了固态电池中一些新的还未用实验验证过的固态电解质。从已知晶体库中选取几千种含锂材料，其中部分材料的固态电解质属性已知，但大部分未知。通过非监督学习方法，鉴定出16个室温离子电导率超过10-4 S cm-1的新型化合物，适合作为固态电解质。其中3个化合物的电导率超过10-2 S cm-1。

近年来，将机器学习与材料学结合是材料科学发展的新趋势。但是，机器学习有相当数量的算法是基于监督学习方法，需要数据来源足够，具有较好的标签。而恰好已知的化合物室温电导率值不完全，且差异较大，想通过监督学习对数千种化合物做出可靠的预测是比较困难的。因此，非监督学习更为合适，因其可以容易地应用于庞大的数据集，而不管是否存在任何属性或标签。

具体方法如下。将各种锂离子化合物阴离子晶格作为特征输入量，获得各种材料的改进XRD图谱（mXRD）。对该图谱进行无监督学习训练，将材料按照电导率性能的优劣进行归类分组。然后对有潜力的分组，再进行精确的第一性原理计算，进一步验证是否是好的固态电解质候选材料。

<br> (二维码自动识别)

首先，在mXRD中，进行无监督学习的clustering(聚类)算法。

首先基于层次聚类法，得到模型C1，发现但大多数RT接近10-3~10-2 S cm-1的固态锂离子导体都聚集在树状图中心的两组中，即组V和组 VI。

此外还依据其他两种不同的聚类技术分别训练了模型C2、C3，结果显示已知的固态锂离子导体在这三个模型生成的组之间有很大的重叠，从而证实了基于mXRD和无监督学习算法能够对快速导电材料与不良导电材料进行聚类分类。

他们这个工作中代码使用的是R语言。

这里我具体介绍下C1模型和代码实现。C1方法是基于层次聚类。源代码如下：

Li-ion-hierarchical-clustering-master

# Hierarchical clustering

mydata=read.csv('Liion_comp_528.csv') 首先载入原始数据。一共380个材料，也就是有380行。纵列中，可以看到材料的名字, ID，电导率cond，电导率的对数，后面V1，V2，V3可以理解成mXRD对应每个角度的强度值。V一共有4504列

这里我以In16Li8S32Sn4为例，将1-4501作为横坐标，V1、V2、V3......V4501作为纵坐标，绘制了X-Y图。显示就是其mXRD图

然后看看mydata读取之后是个什么东西

> attributes(mydata)

$names

[1] "formula_id" "index" "cond" "log10_cond" "V1"

[6] "V2" "V3" "V4" "V5" "V6"…

[996] "V992" "V993" "V994" "V995" "V996"

$class

[1] "data.frame"

X_ini=as.matrix(mydata[,5:4505])

# known li-ion conductivity values

Mydata本来是data.frame格式

读取从第5列到4505列的列数据，并转成矩阵形式X_ini

来看下X_ini的属性

> attributes(X_ini)

$dim

[1] 528 4501

$dimnames

$dimnames[[1]]

NULL

$dimnames[[2]]

[1] "V1" "V2" "V3" "V4"

y=mydata$log10_cond

读取了第四列数据

# Examine number of clusters from hierarchical clustering, by examine the variance# function returning the variance
myobj=function(X){
ncut_vec=seq(2,10)

hc_df=mydata[,c('formula_id','index','log10_cond')]

取3列形成新的数据hc_df

hc_df[,as.character(ncut_vec)]=NA

$names

[1] "formula_id" "index" "log10_cond" "2"

[5] "3" "4" "5" "6"

[9] "7" "8" "9" "10"

$class

[1] "data.frame"

等于是强行在后面加了2-10的colnames，然后内容都是NA

原本有数字的部分，依然还是有的

hc_run <- hclust(dist(X),method="ward.D2")

我用X_ini作为X运行了一次

这里，X_ini是一个矩阵

结果如下

> attributes(hc_run)

$names

[1] "merge" "height" "order" "labels"

[5] "method" "call" "dist.method"

$class

[1] "hclust"

> hc_run

Call

hclust(d = dist(X_ini), method = "ward.D2")

Cluster method : ward.D2

Distance : euclidean

Number of objects: 528

这里先图形化一下

> plot(hc_run)

关于层次聚类法和hclust的用法，这里说明下。

层次聚类法。先计算样本之间的距离。每次将距离最近的点合并到同一个类。然后，再计算类与类之间的距离，将距离最近的类合并为一个大类。不停的合并，直到合成了一个类。其中类与类的距离的计算方法有：最短距离法，最长距离法，中间距离法，类平均法等。比如最短距离法，将类与类的距离定义为类与类之间样本的最段距离。。。

r语言中使用hclust(d, method = "complete", members=NULL) 来进行层次聚类。

其中d为距离矩阵。

method表示类的合并方法，有：

single 最短距离法

complete 最长距离法

median 中间距离法

mcquitty 相似法

average 类平均法

centroid 重心法

ward 离差平方和法

for (cut_dex in 1:length(ncut_vec)){
ncuts=ncut_vec[cut_dex] 取序列seq(2,10)中的整数
mycut=cutree(hc_run,ncuts) 取出的整数会作用于mycut并赋值给hc_df

cutree函数。ncut代表被切割成组的数量，或者切割高度。

按照官方说法，cutree返回的是：

cutree returns a vector with group memberships if k or h are scalar, otherwise a matrix with group memberships is returned where each column corresponds to the elements of k or h, respectively (which are also used as column names).

运行下

> length(ncut_vec)

[1] 9

假定cut_dex我取为2

> ncuts=ncut_vec[2]

> ncuts

[1] 3

> mycut=cutree(hc_run,ncuts)

hc_run本来就是之前已经分好类的，能plot出来图像的树状图结果。这里再被cutree变成 mycut

> mycut

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[30] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2

[59] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

[88] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3

[117] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

......

[436] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

[465] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

[494] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

[523] 3 3 3 3 3 3

尝试将ncuts从3换成7

> mycut2=cutree(hc_run,7)

> mycut2

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[30] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2

[59] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

[88] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 4 4 4 4 4 4 4 4

[117] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4

[146] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5

[175] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

[204] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

[233] 5 5 5 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6

[262] 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

[291] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

......

[494] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7

[523] 7 7 7 7 7 7

hc_df[,as.character(ncuts)]=mycut
} 这里X和hc_df发生相互作用的点在于，X被聚类的结果hc_run会变成mycut再赋值给hc_df

hc_df的每一列都是聚类之后的hc_run按照不同的ncuts值cutree的结果。如第二列就是ncuts=2的结果，第三列就是ncuts=3的结果

具体案例如下

> ncuts

[1] 3

> as.character(ncuts)

[1] "3"

> hc_df[,as.character(3)]

[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

[20] NA NA……

[514] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

> hc_df[,3]

[1] NA NA -6.080922 NA NA

[6] NA -6.236572 -6.795880 NA NA

[11] NA……

[526] NA NA -14.236572

可见hc_df[,3]和hc_df[,”3”]是不同的

前者其实是hc_df$"log10_cond"，后者是hc_df$”3”

hc_df[,as.character(ncuts)]=mycut

左右两个都有528个元素

等于把11112222333…赋值给NANANANANA

之后，继续attributes

> attributes(hc_df)

$row.names

[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14

[15] 15 16 17

$names

[1] "formula_id" "index" "log10_cond" "2"

[5] "3" "4" "5" "6"

[9] "7" "8" "9" "10"

$class

[1] "data.frame"

下面，整体把上面的for循环跑完

> for (cut_dex in 1:length(ncut_vec)){

+ ncuts=ncut_vec[cut_dex]

+ mycut=cutree(hc_run,ncuts)

+ hc_df[,as.character(ncuts)]=mycut

+ }

结果如下

> hc_df

formula_id index log10_cond 2 3 4 5 6 7 8 9 10

1 In16Li8S32Sn4 1 NA 1 1 1 1 1 1 1 1 1

2 Li1O2Ti1 2 NA 1 1 1 1 1 1 1 1 1

3 rocksalt 3 -6.080922 1 1 1 1 1 1 1 1 1

4 Br4Li2Mg1 4 NA 1 1 1 1 1 1 1 1 1

> hc_df[522:528,]

formula_id index log10_cond 2 3 4 5 6 7 8 9 10

522 Ba2Li3O21P7 522 NA 2 3 4 5 6 7 8 9 10

523 C1Ba9Cl7Li1O28Si10 523 NA 2 3 4 5 6 7 8 9 10

524 B1Li1O14S4 524 NA 2 3 4 5 6 7 8 9 10

525 C3F9Li1O9Rb2S3 525 NA 2 3 4 5 6 7 8 9 10

526 B7Li3O12 526 NA 2 3 4 5 6 7 8 9 10

527 B11Li3O18 527 NA 2 3 4 5 6 7 8 9 10

528 B3Li2O8P1 528 -14.23657 2 3 4 5 6 7 8 9 10

看起来mycut是都录入到了hc_df中

hc_sum=data.frame(matrix(NA,nrow=length(ncut_vec),ncol=1))

> length(ncut_vec)

[1] 9

建立9行1列的空矩阵，并转化成 df

> attributes(hc_sum)

$names

[1] "matrix.NA..nrow...length.ncut_vec...ncol...1."

$class

[1] "data.frame"

$row.names

[1] 1 2 3 4 5 6 7 8 9

> hc_sum

matrix.NA..nrow...length.ncut_vec...ncol...1.

1 NA

2 NA

3 NA

4 NA

5 NA

6 NA

7 NA

8 NA

9 NA

colnames(hc_sum)=c('label_var')rownames(hc_sum)=as.character(ncut_vec)# use the measure of summing both labeled and unlabeled

> hc_sum

label_var

2 NA

3 NA

4 NA

5 NA

6 NA

7 NA

8 NA

9 NA

10 NA

for (group_meth in as.character(ncut_vec)){
aggcount=aggregate(hc_df[,'log10_cond'],list(hc_df[,group_meth]),function(x)sum(!is.na(x)))

令group_mech = “3“

> group_meth <- as.character(3)

> group_meth

[1] "3"

> hc_df[,group_meth]

[1] 1 1 1 1 1

[494] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

[523] 3 3 3 3 3 3

可以搜索下aggregate的用法。这里的函数是对!is.na(x)部分进行求和

is.na(x)判断是否缺失

运行结果

>aggcount=aggregate(hc_df[,'log10_cond'],list(hc_df[,group_meth]),function(x)sum(!is.na(x)))

> aggcount

Group.1 x

1 1 14

2 2 13

3 3 54

相当于1大类里面，有14个材料的'log10_cond'非空

out_stat=matrix(NA,nrow=nrow(aggcount),ncol=3)colnames(out_stat)=c('gp','count','var')

> out_stat

gp count var

[1,] NA NA NA

[2,] NA NA NA

[3,] NA NA NA

aggvar=aggregate(na.omit(hc_df)[,'log10_cond'],list(na.omit(hc_df)[,group_meth]),var)

na.omit（<向量a>）: 返回删除NA后的向量a

函数是var，var：协方差阵及相关阵计算

结果

> aggvar

Group.1 x

1 1 3.296540

2 2 9.491358

3 3 17.256080

out_stat[,'gp']=aggcount[,1]
out_stat[,'count']=aggcount[,2]
out_stat[match(aggvar[,1],out_stat[,'gp']),'var']=aggvar[,2]

> aggvar[,1]

[1] 1 2 3

> out_stat[,'gp']

[1] 1 2 3

> match(aggvar[,1],out_stat[,'gp'])

[1] 1 2 3

myout3=sum(out_stat[,'count']*out_stat[,'var'],na.rm=T)

第二列数量和第三列协方差乘积，再求和

> myout3

[1] 1101.368

hc_sum[group_meth,'label_var']=myout3
}

for循环跑完的结果

> hc_sum

label_var

2 1090.2884

3 1101.3675

4 1088.2601

5 1089.3562

6 1099.3131

7 859.2089 确实从6到7的drop最大。如果cluster成6个，不准，cluster成8个，浪费资源和算力。7个正好

8 858.7213

9 856.2518

10 835.9830

return((hc_sum))
}# print the variance for different number of cutsmyobj(X_ini)# cut into 7 clusters, since it gives the largest drop in variance
hc_out=hclust(dist(X_ini),method='ward.D2')
mycut=cutree(hc_out,7)

> mycut=cutree(hc_out,7)

> mycut

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

[30] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2

[59] 2 2 2 2

......

[523] 7 7 7 7 7 7

图2. 固态锂离子导体无监督学习的示意图

C1模型中，层次聚类，分成7个组的结果如上。

C2模型叫Spectral_clustering 谱聚类。是基于PCA，主成分分析和降维的思维的算法。要比C1复杂很多。代码长度也是C1方法的好几倍。这里不再细，感兴趣的可以去看源代码

C3模型是Subset Optimization

通过这三个模型和各自的聚类算法，从ICSD筛选最初的2989个化合物进行筛选分组，得到各自的快速导电材料集，发现他们三个之间有很大重叠。对这三个集取交集，得到82个独特的化合物。这个无监督学习成功减少了对数千种化合物进行高通量筛选的工作。之后再对这82个化合物进行精确的第一性原理计算，进一步验证固态锂离子的导通性。

图3. 新预测的与已知的固态锂离子导体的离子导电性能对比

AIMD模拟（如图3），结果预测有16个候选化合物的室温电导率高于10-4 S cm-1。尽管这些新发现材料的综合性能不如现有的已知固体电解质材料，但是通过计算可以筛选出最有性能的材料，并发现某些材料在某些方面可以有所改进。

总结与展望

该工作通过使用三种无监督学习模型并取集合的方法，成功地在已知晶体库中，区分了锂离子固体快离子导体和导锂较差的材料。发现其中有16种材料的室温电导率高于10-4 S cm-1，甚至其中一些化合物离子电导率超过10-2 S cm-1。这些新发现的候选物与当前已知的快锂离子导体材料具有截然不同的结构和化学组成，证明了无监督学习方法在新材料发现方面的可行性和高效性。

后面还会介绍更多AI，机器学习，深度学习和材料学结合的优质文章和代码。敬请关注。